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Abstract 

Storm  enhanced  densities  (SEDs)  are  ionospheric  plasma  enhancements  that  disrupt 
radio  communications  in  the  near-Earth  space  environment,  degrading  the  Global 
Positioning  System  (GPS)  and  other  key  technologies.  Accurate  GPS/total  electron 
content  (TEC)  correction  maps  produced  by  ionosphere  models  can  mitigate  degra¬ 
dations  from  SEDs.  An  artificial  SED  was  created  and  ingested  via  slant  TEC  mea¬ 
surements  into  the  Global  Assimilation  of  Ionospheric  Measurements  Gauss-Markov 
Kalman  Filter  Model  to  determine  how  many  ground  GPS  receivers  are  needed  to 
produce  reliable  GPS/TEC  correction  maps  over  the  continental  United  States  during 
geomagnetic  storming.  It  was  found  that  110  well-positioned  GPS  receivers  produced 
the  best  overall  TEC  accuracy,  although  significantly  improved  accuracy  was  still 
achieved  if  40  or  more  receivers  were  used.  Furthermore,  receiver  positioning  had  a 
greater  impact  on  TEC  accuracy  than  the  number  of  receivers  used.  It  was  also  found 
that  TEC  accuracy  for  the  SED  region  increased  at  the  expense  of  TEC  accuracy  ev¬ 
erywhere  else  on  the  map. 
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EFFECTS  OF  STORM  ENHANCED  DENSITIES  ON  GEO-LOCATION 


ACCURACY  OVER  CONUS 

I.  Introduction 

1.1  Motivation  and  Background 

In  recent  years  society  has  greatly  benefitted  from  advances  in  spaced-based  com¬ 
munication  and  navigation  technology.  Civilian,  military,  and  government  institu¬ 
tions  alike  have  thrived  on  the  conveniences  offered  by  these  systems  and  have  come 
to  rely  on  them  daily.  The  Global  Positioning  System  (GPS)  is  one  such  modern 
convenience  that  has  brought  enormous  improvements  to  navigation  and  geo-location 
since  its  creation  in  the  late  1970’s. 

As  with  other  platforms  that  dwell  in  the  near-Earth  environment,  GPS  satellites 
broadcast  signals  to  Earth  via  electromagnetic  transmissions  that  must  pass  through 
a  deep  atmospheric  layer  known  as  the  ionosphere.  The  ionosphere  is  the  layer  of 
the  Earth’s  atmosphere  consisting  of  free  thermal  electrons  and  ions,  as  well  as  some 
neutral  species.  It  extends  from  about  50  -  1500  km  and  then  transitions  into  the 
plasmasphere  ( Schunk  and  Nagy ,  2009).  Several  regions  exist  within  the  ionosphere, 
the  characteristics  of  which  are  highly  dependent  on  time  of  day,  season,  solar  cycle, 
geographic  location,  and  geomagnetic  activity. 

The  ionosphere  has  significant  spatial  and  temporal  variations  that  can  enhance, 
degrade,  or  completely  disrupt  electromagnetic  transmissions.  Because  of  society’s 
increased  dependence  on  GPS  and  other  systems  that  employ  electromagnetic  prop¬ 
agation  through  the  ionosphere,  the  need  for  accurate  ionospheric  forecasts  has  also 
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increased.  Reliable  space  weather  forecasts  can  help  mitigate  degradation  to  impor¬ 
tant  military,  government,  and  civil  platforms  caused  by  ionospheric  disturbances. 

Geomagnetic  storming  typically  causes  large  disturbances  in  the  ionosphere  that 
may  last  for  several  days  ( Prolss ,  2004).  Although  the  most  active  periods  of  geomag¬ 
netic  storming  occur  during  the  peak  of  the  11-year  solar  cycle,  storming  can  happen 
at  solar  minimum  as  well  due  to  solar  phenomena  like  coronal  mass  ejections  (CMEs) 
and  corotating  interactive  regions  ( Schunk  and  Nagy ,  2009). 

Regions  of  increased  plasma  density  known  as  storm  enhanced  densities  (SEDs) 
are  regularly  observed  during  geomagnetic  storms  and  can  seriously  degrade  geo¬ 
location  accuracy  by  interfering  with  GPS  signals.  One  system  prone  to  geo-location 
disruption  due  to  SEDs  is  the  Federal  Aviation  Administration’s  Wide  Area  Augmen¬ 
tation  System  (WAAS).  This  air  navigation  aid  augments  GPS  accuracy  and  enables 
precision  approaches  to  any  airport  within  its  coverage  area  ( Federal  Aviation  Admin¬ 
istration,  2011).  For  example,  a  halo  CME  on  3  April  2010  and  ensuing  geomagnetic 
storm  caused  especially  poor  satellite  tracking  at  several  WAAS  sites  in  Alaska,  re¬ 
sulting  in  a  34-minute  Precision  Approach  service  outage  {Rosen,  2010).  Such  outages 
could  prove  dangerous  for  major  airports  during  hours  of  heavy  air  traffic. 

Ionosphere  models  reduce  the  effects  of  geomagnetic  storming  on  geo-location 
accuracy  by  providing  GPS/Total  Electron  Content  (TEC)  correction  maps,  enabling 
GPS  receivers  to  correct  for  ionospheric  enhancements  that  alter  GPS  signals.  The 
model  used  operationally  by  the  Air  Force  Weather  Agency  (AFWA)  is  the  Global 
Assimilation  of  Ionosphere  Measurements  Gauss-Mar kov  (GAIM-GM)  Kalman  Filter 
Model,  developed  by  Dr.  Robert  Schunk  and  his  team  at  Utah  State  University 
(USU).  This  model  uses  the  physics-based  Ionospheric  Forecast  Model  (IFM)  and 
data  assimilation  from  a  variety  of  sources  to  create  accurate  specifications  of  the 
ionosphere.  By  ingesting  real  observations  of  the  ionosphere,  GAIM-GM  has  been 
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shown  to  create  TEC  maps  that  significantly  improve  the  comparison  to  independent 
TEC  observations  [( Thompson  et  al. ,  2005)  ( Scherliess  et  al. ,  2006),  ( Zhu  et  al. ,  2006), 
(Decker  and  McNamara ,  2007),  (Schunk  et  al ,  2010)]. 

1.2  Research  Objective 

The  objective  of  this  research  is  to  determine  how  many  ground  GPS  receivers  (i.e. 
ground  stations)  are  needed  in  the  continental  United  States  (CONUS)  in  order  to 
have  reliable  GPS/TEC  correction  maps  during  geomagnetic  storming.  This  is  done 
by  imposing  an  artificial  SED  on  an  electron  density  background  and  then  taking  slant 
TEC  measurements  through  that  background  via  simulated  GPS  satellites  and  several 
different  distributions  of  GPS  ground  stations.  GAIM-GM  ingests  these  data  and 
produces  TEC  maps  that  differ  in  accuracy  according  to  the  number  of  GPS  ground 
stations  used.  By  varying  the  number  of  available  ground  stations  and  comparing 
the  model  output,  a  clearer  picture  emerges  of  how  many  ground  GPS  receivers  are 
needed  for  reliable  TEC  maps. 

1.3  Document  Structure 

This  document  is  organized  into  five  chapters.  Chapter  II  provides  background  in¬ 
formation  on  the  ionosphere  and  SEDs  caused  by  geomagnetic  storming.  It  also  gives 
an  overview  of  the  GAIM-GM  model.  Chapter  III  explains  the  methodology  used  to 
conduct  this  research  and  outlines  the  selection  of  model  parameters,  designing  the 
synthetic  SED,  piercing  the  SED  with  slant  TEC  measurements,  running  GAIM-GM, 
and  data  comparison  techniques.  Chapter  IV  presents  the  results  of  assimilating  slant 
TEC  data  into  GAIM-GM  for  varying  GPS  ground  station  distributions  and  gives  skill 
score  analyses.  A  summary  of  the  results,  final  conclusions,  and  recommendations  for 
future  work  are  given  lastly  in  Chapter  V. 
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II.  Background 


2.1  Chapter  Overview 

This  chapter  provides  background  information  on  the  ionosphere  and  geomagnetic 
storms,  and  details  how  these  storms  impact  the  GPS.  The  basic  structure  of  the 
ionosphere  and  the  processes  by  which  it  is  formed  are  described,  followed  by  an 
overview  of  ionospheric  disturbances  known  as  SEDs  and  their  effects  on  human 
activity.  The  two  ionosphere  models  used  in  this  research  are  also  described,  namely, 
the  IFM  and  the  GAIM-GM  Kalman  Filter  Model. 

2.2  The  Ionosphere 

The  ionosphere  is  the  region  of  the  Earth’s  atmosphere  where  there  are  enough 
free  thermal  electrons  and  ions  to  create  a  plasma  dense  enough  to  affect  radio  wave 
propagation.  Some  of  these  ions  and  free  electrons  are  produced  by  high  latitude 
particle  precipitation  but  the  majority  result  from  photoionization  of  neutral  species 
in  the  upper  atmosphere  (above  50  km)  by  solar  and  stellar  irradiance  in  the  extreme 
ultraviolet  (EUV)  and  x-ray  wavelength  range. 

The  photoionization  rate  in  Earth’s  atmosphere  can  be  represented  by  the  highly 
simplified  Chapman  production  function: 

P{z,x)  =  I{z,x)rivan(z)  (1) 

In  this  equation  P  is  the  ion  production  rate ,  /  is  the  photon  flux  as  a  function  of 
height  z  and  solar  zenith  angle  y,  rj  is  the  probability  of  photon  absorption  (resulting 
in  the  production  of  an  electron-ion  pair),  aa  is  the  absorption  cross-section ,  and 
n(z)  is  the  neutral  species  density  as  a  function  of  height  ( Schunk  and  Nagy ,  2009). 
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Figure  1  illustrates  the  relationship  between  neutral  density  and  solar  irradiance  in 
determining  the  height  of  maximum  photoionization,  known  as  a  Chapman  layer. 


Figure  1.  Chapman  layer  and  ionization  curve  (dashed  line).  Decreasing  neutral  den¬ 
sities  (n(z))  and  increasing  solar  irradiance  (I(z))  with  altitude  form  an  ionization  peak 
known  as  a  Chapman  layer 


As  plasma  is  produced  it  gets  transported  throughout  the  atmosphere  via  different 
mechanisms  such  as  electromagnetic  forces,  gravitational  force,  the  pressure  gradient, 
and  collisions  with  neutral  particles  ( Ondoh  and  Marubashi ,  2001).  Combined  with 
ion/neutral  production  and  loss  processes  through  chemical  reactions,  several  distinct 
layers  of  the  ionosphere  are  formed:  1)  the  D  region;  2)  the  E  region;  3)  the  F  region; 
and  4)  the  topside  ionosphere,  which  transitions  into  the  plasmasphere  at  about  1500 
km.  The  vertical  profile  of  these  regions  are  illustrated  in  Figure  2,  which  also  includes 
typical  ion  composition  and  electron  density  values  for  each  region. 

The  lowest  layer  of  the  ionosphere,  the  D  region,  is  found  between  50  -  90  km. 
It  is  the  most  complex  of  all  four  regions  and  here  ion  production  and  loss  are  dom¬ 
inated  by  chemical  processes  occurring  in  the  relatively  dense  neutral  atmosphere. 
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Figure  2.  Vertical  profile  of  the  ionosphere.  Densities  consisting  mainly  of  0+  reach 
maximum  values  near  the  F2  peak,  seen  at  about  300  km  in  this  figure.  Additionally, 
a  smaller  peak  consisting  of  several  ion  species  is  found  in  the  E  region  near  130  km 
(Adapted  from  S chunk  and  Nagy  (2009)) 


Molecular  ions  are  most  abundant  in  this  region  and  there  also  exist  negative  ions 
which  contribute  significantly  to  total  plasma  density.  Above  this  layer,  the  E  region 
spans  from  about  90  -  150  km.  Its  dominant  ions,  NO+  and  Oj,  are  formed  from 
direct  ionization  of  02,  as  well  as  from  chemical  processes.  High  electron  densities  in 
this  region  do  not  occur  because  of  an  increased  rate  of  ion  recombination  reactions. 

Above  the  E  region  is  the  F  region,  which  consists  of  three  sub-regions:  1)  Fi 
region;  2)  F2  region;  and  3)  upper  F  region  or  topside  ionosphere.  Transport  processes 
in  the  F  region,  some  of  which  are  very  complex,  are  more  dominant  than  those  in  the 
D  and  E  regions  and  govern  much  of  the  ion  production  and  loss.  The  Fi  region  is 
found  between  approximately  150  -  250  km.  Plasma  transport  in  the  Fx  region  leads 
to  higher  electron  densities  here  than  in  lower  regions  and  these  densities  continue  to 
increase  until  peaking  in  the  F2  region.  This  region  spans  from  about  250  -  600  km 
and  contains  the  highest  electron  densities  of  the  ionosphere.  The  F2  peak  can  be 
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seen  in  Figure  2  near  the  300  km  height  where  electron  densities  approach  106  cm~3. 
Above  the  F2  region,  the  layer  extending  from  about  600  -  1500  km  is  referred  to  as 
the  topside  ionosphere.  The  region  of  dominant  H+  ions  above  the  topside  ionosphere 
is  called  the  plasmasphere  (or  protonosphere)  (S chunk  and  Nagy ,  2009). 

The  ionosphere  is  often  measured  in  terms  of  TEC  where  one  TEC  unit  equals 
1016  electrons / m2 .  It  represents  the  amount  of  free  electrons  present  along  a  path 
between  two  points  (typically  the  surface  of  the  Earth  and  the  top  of  the  ionosphere). 
TEC  analyses  can  easily  be  overlaid  onto  world  or  regional  map  projections,  providing 
a  fast  way  to  visualize  the  ionosphere’s  density.  A  more  detailed  explanation  of  how 
TEC  can  be  derived  using  GPS  satellites  is  given  in  Appendix  1. 

All  regions  of  the  ionosphere  are  subject  to  variation  as  a  function  of  latitude,  lon¬ 
gitude,  universal  time,  season,  solar  cycle,  and  geomagnetic  activity.  As  a  result,  the 
ionosphere  is  typically  considered  to  have  two  states:  a  background  or  climatological 
state  and  a  disturbed  or  storm  state  (S chunk  and  Scherliess ,  2004).  Both  states  can 
affect  electromagnetic  wave  propagation,  although  the  disturbed  state  tends  to  cause 
greater  disruptions  due  to  the  effect  of  geomagnetic  storming  on  plasma  densities. 

2.3  Geomagnetic  Storms 

Geomagnetic  storms  are  events  of  strongly  enhanced  dissipation  of  solar  wind 
energy  into  the  near- Earth  space  environment.  During  these  storms,  the  rate  of 
energy  transfer  to  the  magnetosphere  can  be  many  times  greater  than  the  quiet 
condition  energy  transfer  rate  from  the  solar  wind.  The  effects  of  this  impulsive 
energy  transfer  -  which  usually  lasts  for  1  -  3  days  -  disturb  the  magnetosphere, 
thermosphere,  ionosphere,  auroral  ovals,  and  the  numerous  current  systems  which 
couple  these  regions  together  ( Prolss ,  2004).  Most  periods  of  geomagnetic  storming 
occur  during  solar  maximum  (i.e.,  the  peak  of  the  11-year  solar  cycle),  although 
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storms  can  happen  during  solar  minimum  as  well  (S  chunk  and  Nagy ,  2009). 

Solar  phenomena  typically  responsible  for  geomagnetic  storming  are  CMEs  and 
corotating  interactive  regions.  A  CME  is  the  expulsion  of  a  large  amount  (1012  - 
1013  kg)  of  gaseous  solar  mass  into  interplanetary  space.  These  clouds  of  plasma  are 
associated  with  strong  magnetic  fields  than  can  be  either  northward  (Bz  positive) 
or  southward-oriented  (Bz  negative).  If  the  held  is  negative,  reconnection  with 
the  Earth’s  northward-oriented  magnetic  held  allows  solar  wind  particles  to  more 
easily  penetrate  into  the  inner  magnetosphere.  CMEs  with  Bz  negative  magnetic 
helds  that  reach  Earth  have  been  frequently  associated  with  the  onset  of  geomagnetic 
storms.  In  the  case  of  corotating  interactive  regions,  fast  solar  wind  from  a  coronal 
hole  interacts  with  the  slower,  surrounding  solar  wind  and  leads  to  compression  of 
the  interplanetary  magnetic  held  (IMF).  If  the  ambient  IMF  already  possesses  a  B^ 
negative  component,  in  can  be  intensihed  to  the  point  where  a  geomagnetic  storm 
occurs  as  the  region  passes  the  Earth  ( Prolss ,  2004). 

There  are  generally  three  phases  in  a  geomagnetic  storm:  1)  initial;  2)  main;  and 
3)  recovery.  The  initial  phase  begins  with  magnetospheric  compression  caused  by  a 
shock  in  the  IMF  (from  a  CME,  etc.).  During  the  main  phase,  energetic  plasma  from 
the  magnetotail  is  injected  into  the  inner  magnetosphere,  intensifying  the  ring  current 
and  inducing  a  magnetic  held  opposed  to  the  Earth’s  own  held  at  the  surface  -  most 
notably  at  equatorial  latitudes.  This  causes  the  Earth’s  magnetic  held  to  weaken  at 
the  surface.  The  recovery  phase  entails  a  gradual  return  of  the  Earth’s  magnetic  held 
to  pre- storm  levels.  This  can  last  up  to  several  days  while  the  enhanced  ring  current 
subsides  and  excess  particles  are  lost  through  various  mechanisms.  Large  storms 
can  signihcantly  modify  the  composition,  density,  and  circulation  of  the  ionosphere- 
plasmasphere  system  on  a  global  scale.  These  changes  can  persist  for  several  more 
days  after  the  geomagnetic  storming  ceases  ( Schunk  and  Nagy ,  2009). 


A  variety  of  geomagnetic  indices  have  been  devised  to  quantify  variations  in  the 
Earth’s  magnetic  held  and  current  systems  during  geomagnetic  storms.  Some  com¬ 
monly  used  indices  are  the  Kp ,  ap,  and  Dst  indices.  A  typical  fluctuation  profile  for 
some  of  these  indices  during  a  major  geomagnetic  storm  is  shown  in  Figure  3. 


Figure  3.  Profile  of  ap  and  Dst  indices  during  the  10  November  1986  geomagnetic  storm. 
The  linearly  scaled,  3-hourly  ap  index  is  simply  a  derivation  of  the  quasi-logarithmic- 
scaled  Kp  index 

The  Dst  index  is  given  in  nano  teslas  (nT)  and  indicates  disturbances  to  the 
Earth’s  magnetic  held  due  to  both  ring  current  enhancements  and  increased  solar  wind 
pressure  on  the  magnetosphere.  The  Kp  and  ap  indices  are  planetary  measurements 
representing  the  strength  the  Earth’s  magnetic  held  at  the  surface,  relative  to  an 
assumed  quiet-day  curve.  The  Kp  index  uses  a  quasi-logarithmic  scale  while  the  ap 
index  is  derived  from  the  Kp  index  and  uses  a  linear  scale  [( Perrone  and  De  Franceschi , 
1998),  ( National  Geophysical  Data  Center ,  2011)].  In  Figure  3  both  ap  and  Dst 
indices  increased  during  the  initial  storm  phase,  followed  by  a  sharper  increase  in  ap 
and  a  simultaneous  drop  in  Dst  during  the  main  phase.  The  recovery  phase  for  this 
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storm  lasted  about  two  days  as  energetic  particles  which  enhanced  the  ring  current 
precipitated  out  of  the  magnetosphere  and  both  indices  returned  to  normal  values. 


2.4  Storm  Enhanced  Density  (SED) 

During  the  early  stages  of  geomagnetic  storms,  regions  of  increased  plasma  density 
are  regularly  observed  equatorward  of  the  auroral  oval  in  the  pre-midnight  ionosphere. 
These  SEDs  consist  of  ionization  plumes  of  enhanced  F  region  plasma  where  TEC  can 
be  two  to  four  times  greater  than  background  TEC  values  ( Foster ,  1993).  SEDs  are 
prevalent  over  North  America  during  high  geomagnetic  activity.  A  list  of  the  defining 
SED  characteristics  and  their  values  are  presented  in  Table  1. 


Table  1.  SED  properties 


Property 

Typical  value 

Source 

Duration 

2-5  hours 

( Coster  et  al. ,  2001),  ( Foster  et  al. ,  2004) 

Peak  time 

1730  -  2000  UT 

(. Foster ,  1993),  ( Coster  et  al. ,  2001),  ( Foster 
et  al. ,  2005),  ( Foster  et  al. ,  2006) 

HmF2 

>  500  km 

(. Foster  et  al. ,  2005) 

Altitude  extent 

300  -  500  km 

(Foster,  1993),  (Foster  et  al. ,  2004) 

Width 

600  -  1000  km 

(Foster  et  al. ,  2004),  (Kelley  et  al. ,  2004) 

Velocity 

0.5  -  1  km/s 

(Foster,  1993),  (Coster  et  al,  2001),  (Foster 
et  al,  2002),  (Foster  et  al,  2005) 

NmF2 

1-3  x  1012  m~J 

(Foster  et  al,  2004),  (Foster  et  al,  2005) 

TEC  increase 

2-4  times 

(Foster,  1993),  (Coster  et  al,  2001) 

Figure  4  illustrates  the  storm-time  convection  that  sets  up  to  form  SEDs.  The 
plasma  found  within  SEDs  actually  originates  in  the  sunlit  equatorial  ionosphere. 
Early  in  the  geomagnetic  storm  a  strong,  eastward  electric  held  at  low  latitudes 
drives  plasma  upward  into  the  plasmasphere  via  an  E  x  B  drift  so  fast  that  it  cannot 
dissociatively  recombine  with  free  electrons.  From  high  altitudes  it  diffuses  down 
geomagnetic  field  fines  to  lower-mid  latitudes.  At  this  point  a  poleward-oriented 
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Eastward  electric  field  uplifts  equatorial  ionosphere 


Poleward  SAPS  electric  field  strips 
away  outer  layers  of  plasmasphere 


Figure  4.  Plasma  convection  pattern  during  geomagnetic  storming  that  leads  to  SED 
formation.  Enhanced  electric  fields  drive  equatorial  plasma  (light  grey  patches)  high 
into  the  plasmasphere,  where  it  diffuses  down  magnetic  field  lines.  At  lower  mid¬ 
latitudes  an  enhanced  polarizing  electric  field  then  erodes  the  plasma  from  the  plasma¬ 
sphere  and  transports  it  poleward  (dark  grey  patches)  in  the  Sub- Auroral  Polarization 
Stream  (SAPS).  (Adapted  from  Schunk  and  Nagy  (2009)) 


electric  field  convects  the  plasma  both  downward  in  altitude  and  towards  the  poles. 
This  electric  held  arises  during  disturbed  conditions  and  exists  to  close  the  circuit  in 
the  region  of  low  conductivity  between  held  aligned  currents  in  the  nighttime  sub- 
auroral  zone.  Across  this  zone  a  region  of  strong  plasma  drift  known  as  the  sub-auroral 
polarization  stream  (SAPS)  is  formed  and  SEDs  in  the  stream  drift  in  a  northwest 
direction  across  the  U.  S.  ( Foster  et  al. ,  2002). 

Figure  5  shows  the  TEC  characteristics  of  a  SED  plume  over  North  America  during 
a  severe  geomagnetic  storm  on  31  March  2001.  This  plot  was  created  using  ionosphere 
TEC  values  collected  from  more  than  120  GPS  ground  stations  ( Foster  et  al. ,  2002). 
In  typical  cases  such  as  this,  high  conductivities  on  the  dayside  prevent  the  poleward 
held  from  extending  deep  into  the  sunlit  sector  and  how  stagnates,  causing  plasma 
to  build  up  in  a  narrow  channel  along  the  dusk  terminator  ( Kelley  et  al. ,  2004). 

SEDs  and  strong  geomagnetic  activity  in  general  cause  a  wide  variety  of  societal 
and  economic  impacts.  Much  of  the  technical  infrastructure,  space-based  assets,  and 
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Figure  5.  SED  plume  over  the  northern  U.S.  at  1930  UT  during  the  31  March  2001 
geomagnetic  storm.  The  50  TEC  unit  contour  is  outlined  in  red  and  shows  the  position 
of  the  SED.  Data  from  over  120  GPS  ground  stations  were  used  to  reconstruct  this 
TEC  map.  (Adapted  from  Foster  (2002)) 


enabling  technology  of  the  U.S.  and  world  is  frequently  vulnerable  to  degradation  and 
occasionally  to  damage  during  geomagnetic  storming  ( Baker  et  al. ,  2008).  Impacted 
areas  include,  but  are  not  limited  to:  1)  navigation  and  surveying  using  the  GPS 
system;  2)  high-frequency  communications;  3)  satellite  tracking;  4)  satellite  lifetime; 
5)  surveillance;  6)  pipeline  integrity;  7)  electrical  power  grids;  and  8)  over-the-horizon 
radar.  Taking  into  consideration  ground  induced  currents  on  U.S.  power  grids  alone,  it 
is  estimated  that  a  major  blackout  similar  in  size  to  the  14  August  2003  northeastern 
U.S. /Canada  blackout  could  cost  $4  to  $10  billion  in  economic  loss  {U.S. -Canada 
Power  System  Outage  Task  Force ,  2004). 

Ionospheric  disturbances  due  to  geomagnetic  storming  pose  a  particular  threat  to 
military  objectives  on  both  the  strategic  and  tactical  level.  For  example,  capabilities 
such  as  satellite  uplink  for  forward- deployed  units  and  GPS-guided  munitions  have 
improved  the  effectiveness  of  military  operations  in  recent  years.  As  decision  makers 
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rely  more  on  technological  and  information  superiority,  the  risk  of  degradation  to 
key  navigation,  communication,  and  weapon  systems  during  geomagnetic  storming 
becomes  greater.  To  mitigate  these  risks,  advanced  computer  models  that  can  accu¬ 
rately  specify  and  forecast  the  ionosphere  have  been  developed  and  will  continue  to 
play  a  vital  role  in  the  civilian,  government,  and  military  arena. 

2.5  GAIM-GM  Kalman  Filter  Model 

The  ionosphere  model  of  interest  for  this  research  is  the  GAIM-GM  Kalman  Filter 
Model,  developed  by  Dr.  Robert  Schunk  and  his  team  at  USU  and  used  operationally 
at  AFWA.  It  may  also  be  run  on-demand  for  research  studies  through  the  NASA 
Community  Coordinated  Modeling  Center.  GAIM-GM  creates  a  three-dimensional, 
time-dependent  evolution  of  the  global  ionosphere  from  90  -  1400  km.  In  recent 
years  the  accuracy  of  the  model  has  been  validated  against  a  variety  of  independent 
observations  for  varying  geophysical  conditions  [( Thompson  et  al. ,  2005)  ( Scherliess 
et  al. ,  2006),  ( Zhu  et  al. ,  2006),  ( Decker  and  McNamara ,  2007),  ( Schunk  et  al. ,  2010)]. 

GAIM-GM’s  strength  as  an  ionospheric  model  comes  from  its  ability  to  assimilate 
ionospheric  measurements  from  several  sources,  which  include:  1)  bottomside  electron 
density  profiles  from  a  variable  number  of  ionosondes;  2)  in  situ  electron  density  from 
four  Defense  Meteorological  Satellite  Program  (DMSP)  satellites;  3)  nighttime  line- 
of-sight  ultraviolet  radiances  measured  by  satellites;  and  4)  slant  TEC  from  a  variable 
number  of  GPS  satellites  ( Scherliess  et  al. ,  2006).  A  schematic  diagram  of  these  data 
sources  are  shown  in  Figure  6. 

GAIM-GM  is  based  on  the  global,  physics-based  Ionosphere  Forecast  Model  (IFM) 
and  a  Gauss-Markov  Kalman  filter  data  assimilation  algorithm.  The  IFM  is  used  as  an 
initial  background  from  which  ionospheric  densities  are  taken,  and  then  perturbations 
are  superimposed  on  that  background  by  the  Kalman  filter,  based  on  available  data 
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Figure  6.  Portrayal  of  some  of  the  data  sources  that  are  available  for  assimilation  into 
GAIM-GM.  (Adapted  from  Schunk  (2004)) 

and  their  errors.  A  brief  overview  of  the  IFM  and  how  it  operates  is  given  in  the  next 
section. 

2.5.1  Ionospheric  Forecast  Model  (IFM). 

The  IFM  is  a  numerical,  physics-based  model  that  provides  the  three-dimensional, 
time- dependent  evolution  of  the  global  ionosphere  from  90  -  1600  km.  It  provides  the 
density  distributions  of  four  major  ions  (NO+,  Oj,  Nj,  0+)  at  E  region  altitudes, 
two  major  (NO+,  0+)  and  two  minor  (Nj,  Oj)  ions  at  F  region  altitudes,  TEC, 
and  also  ion  and  electron  temperatures  at  both  E  and  F  region  altitudes.  The  IFM 
also  calculates  H+  densities  in  the  F  region  and  topside  ionosphere,  as  well  as  the 
density  of  the  E  and  F  region  peaks  and  their  altitudes.  These  peaks  are  easily 
resolved  by  the  model’s  vertical  resolution  of  4  km  in  the  E  region  and  20  km  in  the  F 
region  ( Schunk  et  al .,  1997).  The  latitude  and  longitude  resolutions  are  3°  and  7.5°, 
respectively,  and  the  output  can  be  in  either  geographic  or  geomagnetic  coordinates 
(Schunk  and  Scherliess ,  2004).  Figure  7  gives  an  example  of  an  IFM  specification 


14 


(a)  TEC 


-80  -60  -40  -20  0  20 

Latitude 


4.5  5  5.5 

Log  Electron  Density 


(b)  Electron  density  longitude  slice,  240°  E 


Figure  7.  IFM  specification  during  solar  minimum  (Year  2007,  Day  355,  2200  UT). 
Note  how  the  model  output  captures  the  equatorial  Appleton  Anomaly,  or  fountain 
effect,  in  the  240°  E  longitude  slice  (b) 


for  the  winter  solstice  during  solar  minimum.  The  TEC  at  2200  UT  is  shown  in  7(a) 
while  an  electron  density  vertical  profile  at  240°  E  longitude  for  the  same  time  is 
shown  in  7(b).  Note  how  the  IFM  captures  the  equatorial  Appleton  Anomaly. 

The  IFM  uses  the  FI 0.7,  Kp ,  and  Ap  indices  as  solar  and  geomagnetic  inputs. 
These  values  are  in  turn  used  to  compute  neutral  densities,  winds,  temperatures,  and 
the  plasma  convection  and  precipitation  patterns  via  integrated  empirical  models 
(Space  Environment  Corporation ,  2002). 

The  most  important  chemical  and  physical  processes  of  the  ionosphere  are  ac¬ 
counted  for  in  the  IFM’s  numerical  solution.  This  solution  is  based  on:  1)  field-aligned 
plasma  diffusion;  2)  cross-held  electro  dynamic  drifts  governed  by  magnetospheric  and 
dynamo  electric  fields;  3)  thermospheric  winds;  4)  exchange  fluxes  in  the  protono¬ 
sphere;  5)  energy-dependent  chemical  reactions;  6)  electron  thermal  conduction;  7) 
numerous  ion  production  sources  (photoionization  from  solar  EUV  radiation,  auroral 
particle  precipitation,  backscattered  solar  radiation,  and  starlight);  8)  neutral  com¬ 
position  changes;  and  9)  a  variety  of  local  heating  and  cooling  processes  (S chunk  and 
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Scherliess ,  2004).  Validations  of  the  IFM  against  observations  for  different  parame¬ 
ters  have  been  performed  by  Scherliess  et.  al.  (2006),  Zhu  et.  al.  (2006),  Decker  and 
McNamara  (2007),  Thompson  et.  al.  (2006),  and  others. 

2.5.2  GAIM-GM  Output  and  Resolution. 

The  output  from  GAIM-GM  is  a  continuous,  time-varying,  three-dimensional  re¬ 
construction  of  the  electron  density  distribution  from  92  -  1400  km.  The  model 
also  produces  auxiliary  parameters  including  the  heights  and  densities  of  the  F2  re¬ 
gion  peak,  the  heights  and  densities  of  the  E  region  peak,  and  TEC.  GAIM-GM 
also  provides  global  distributions  for  the  ionospheric  drivers,  which  include  magne- 
tospheric  and  equatorial  electric  fields,  electron  precipitation  patterns,  and  neutral 
winds  ( Schunk  and  Scherliess ,  2004). 

There  are  two  modes  for  running  the  model:  global  and  regional.  Figure  8  shows 
the  latitudinal  and  longitudinal  resolutions  of  both  modes.  The  global  mode  has  a 
typical  resolution  of  4.6°  latitude  by  15°  longitude,  while  the  regional  mode  is  user- 
adjustable  and  can  resolve  up  to  1°  latitude  by  3.75°  longitude,  based  on  data  coverage 
and  the  computational  environment.  The  vertical  resolution  of  the  output  is  4  km 
in  the  E  region  and  20  km  in  the  F  region  and  above,  for  both  global  and  regional 
modes  ( Scherliess  et  al. ,  2006). 

2.5.3  Gauss-Mar kov  Kalman  Filter. 

GAIM-GM  operates  by  using  a  Kalman  filter  to  superimpose  density  perturbations 
onto  background  IFM  density  field,  based  on  available  data  and  their  errors.  These 
perturbations  and  errors  then  evolve  over  time  with  a  statistical  model  (the  Gauss- 
Markov  process).  The  total  electron  density  at  each  grid  point  is  expressed  as: 

N  =  NIFM  +  Npert  (2) 
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Figure  8.  GAIM-GM  model  output  for  global  and  regional  modes.  Global  mode 
resolution  is  typically  4.6°  latitude  by  15°  longitude,  while  resolution  as  high  as  1° 
latitude  by  3.75°  longitude  can  be  obtained  in  regional  mode 


where  NIFM  is  the  electron  density  from  the  IFM  background  and  Npert  is  the  pertur¬ 
bation  density  derived  by  the  Kalman  filter.  The  density  Npert  is  given  in  a  geographic 
frame  and  evolves  over  each  15-minute  assimilation  time  step  through  the  following 
formula: 

Ktrt  =  L  Nlrt  (3) 

where  L  is  a  transition  matrix,  the  product  of  a  translation  matrix  Li  and  a  diagonal 
matrix  L2.  The  matrix  Li  moves  the  perturbation  density  held  at  each  time  step  in  a 
sun-synchronous  magnetic  frame  while  the  diagonal  matrix  L2  relaxes  the  perturba¬ 
tions  to  zero  if  there  is  no  data  present.  The  relaxation  time  r  is  set  to  r  =  5  hours 
by  default,  but  can  be  adjusted  to  better  capture  changing  geophysical  conditions 
(. Scherliess  et  al. ,  2006). 

The  Gauss-Markov  Kalman  filter  uses  an  error  covariance  matrix  that  evolves 
with  the  same  transition  matrix  L  as  the  electron  density  perturbations.  Included  in 
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the  error  covariance  matrix  are  the  uncertainties  of  the  L,  which  were  obtained  from 
1107  individual  2-day  runs  of  the  IFM  with  varying  external  forcing  parameters  and 
a  wide  range  of  climatological  variations.  In  the  absence  of  assimilation  data  (GPS, 
DMSP,  etc.)  the  model  error  covariances  are  given  by  these  uncertainties.  For  this 
reason,  GAIM-GM  will  not  exactly  reproduce  the  IFM  density  background  even  if  no 
observational  data  are  assimilated. 

Quality  control  is  performed  on  the  observations  assimilated  by  the  Kalman  filter. 
Obviously  wrong  observations  are  rejected  and  appropriate  data  errors  are  assigned 
to  remaining  data.  These  errors  consist  in  an  observational  error  associated  with  the 
instrument  taking  the  measurements,  and  an  error  associated  with  the  representative¬ 
ness  of  the  observation.  All  data  errors  are  assumed  to  have  a  Gaussian  distribution 
and  to  be  uncorrelated  with  one  another  ( Scherliess  et  al ,  2006).  A  problem  can  arise 
when  closely  spaced  measurements  are  assimilated,  such  as  those  from  a  dense  grid  of 
GPS  ground  stations.  The  error  of  representativeness  becomes  highly  correlated  and 
the  TEC  measurements  are  trusted  too  much.  The  solution  available  to  the  user  is  to 
reduce  the  data  density,  either  by:  1)  increasing  the  model  grid  density  via  running 
GAIM-GM  in  regional  mode  with  higher  resolution;  or  2)  maintaining  the  same  reso¬ 
lution  but  simply  using  less  data.  Increasing  the  number  of  ground  stations  without 
raising  grid  resolution  does  not  improve  the  TEC  result  and  can  actually  degrade  the 
F2  region  peak  density  results  ( Thompson  et  al. ,  2005). 

2.5.4  Assimilation  of  GPS  Slant  TEC. 

Among  the  many  data  types  that  GAIM-GM  can  ingest,  GPS  slant  TEC  mea¬ 
surements  are  some  of  the  most  readily  available.  For  a  15-minute  sampling  cadence, 
a  single  GPS  ground  station  located  within  CONUS  can  produce  between  35  -  45 
slant  TEC  measurements  per  hour,  yielding  about  1000  measurements  per  day.  This 
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becomes  a  truly  vast  reservoir  of  data  for  GAIM-GM  when  multiplied  by  the  number 
of  available  GPS  ground  stations. 

The  Gauss-Markov  Kalman  filter  makes  an  important  adjustment  to  slant  TEC 
data  during  the  assimilation  process.  Recall  that  slant  TEC  measurements  consist  of 
integrated  electron  densities  along  the  paths  from  GPS  satellites  to  ground  stations. 
These  satellites  orbit  near  20,200  km  and  thus  slant  TEC  paths  include  electron  con¬ 
tributions  from  both  the  ionosphere  and  plasmasphere  (see  Figure  9).  When  satellites 
are  not  directly  overhead,  slant  paths  can  be  considerably  larger,  causing  a  significant 
amount  of  plasmasphere  electrons  to  be  integrated  into  the  total  slant  path.  Since 
the  upper  bounds  of  GAIM-GM  extend  no  higher  than  1400  km ,  the  plasmaspheric 
contribution  to  GPS  slant  TEC  measurements  (1400  -  20,200  km)  must  be  accounted 
for  and  discarded.  This  is  done  through  an  internal  representation  of  the  plasmas¬ 
phere  which  estimates  the  plasmasphere  TEC  contribution  and  subtracts  it  during 
the  data  assimilation  process.  Even  this  simple  representation  of  the  plasmasphere 
provides  significantly  better  distributions  of  F2  peak  density  and  TEC,  compared  to 
results  where  the  plasmasphere  contribution  is  ignored  ( Thompson  et  al. ,  2009). 

After  the  slant  TEC  data  have  been  adjusted  for  plasmasphere  contributions,  the 
Kalman  filter  assigns  errors  to  the  observations.  One  of  these  errors  is  the  phase 
TEC  error.  This  results  from  inherent  fluctuations  in  the  slant  TEC  estimate  due 
to  bandwidth  limitation  precision,  multipath  errors,  receiver  noise,  and  other  effects 
(. Dyrud  et  al. ,  2008).  Fortunately  the  phase  TEC  error  values  are  normally  given 
separately  in  slant  TEC  observations,  enabling  the  Kalman  filter  to  appropriately 
weigh  ingested  data  when  constructing  the  perturbation  density  matrix  Npert. 

Additionally,  the  Kalman  filter  assigns  errors  according  to  instrument  biases  in 
both  GPS  satellite  and  receiver  sensors.  These  so-called  differential  code  biases 
(DCBs)  are  caused  by  delay  and  dispersion  of  GPS  signals  traveling  through  ana- 
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Figure  9.  Plasmaspheric  contribution  to  slant  TEC  path  between  a  GPS  ground  station 
and  GPS  satellite.  The  GAIM-GM  Kalman  filter  uses  an  internal  algorithm  to  discard 
all  extra  slant  TEC  above  the  data  assimilation  region  of  1400  km 

log  circuits  and  filters  in  satellites  and  receivers,  and  therefore  DCBs  are  unique  to 
each  piece  of  equipment.  DCBs  cause  significant  errors  in  TEC  estimates  and  must  be 
addressed  before  the  data  can  be  reliably  used  for  ionospheric  reconstruction.  Fortu¬ 
nately,  the  International  Global  Navigational  Satellite  System  Service  (IGS)  compiles 
a  daily  list  of  corresponding  satellite-to-receiver  DCBs  for  many  stations  worldwide 
( Dyrud  et  al. ,  2008).  GAIM-GM  ingests  these  lists  from  the  Internet  and  provides  the 
DCBs  to  the  Kalman  filter.  If  GPS  DCBs  are  not  available  for  a  given  ground  station 
the  Kalman  filter  solves  internally  for  the  missing  bias  ( Scherliess  et  al. ,  2006). 

By  assimilating  slant  TEC  from  available  of  ground  GPS  receivers,  GAIM-GM  can 
significantly  improve  its  ionosphere  specification  compared  to  independent  measures 
of  vertical  TEC  and  thus  create  more  reliable  ionosphere  reconstructions  (S chunk 
et  al. ,  2010).  This  will  lead  to  better  ionosphere  correction  maps  and  reduced  degra¬ 
dation  to  military,  government,  and  civil  navigation  and  communication  systems. 
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III.  Methodology 


3.1  Chapter  Overview 

This  chapter  describes  the  methods  used  to  perform  this  study  and  is  broken 
into  seven  sections.  First  it  details  how  the  IFM  density  backgrounds  were  created, 
followed  by  a  section  that  describes  how  the  synthetic  SED  was  modeled  and  imposed 
onto  the  IFM  background.  The  GPS  ground  station  distributions  and  the  methods 
for  creating  them  will  be  examined,  followed  by  a  description  of  how  slant  TEC 
measurements  were  cut  through  the  IFM  background  and  prepared  for  assimilation 
into  GAIM-GM.  This  leads  directly  into  an  outline  of  procedures  for  running  the 
GAIM-GM  model.  Additionally,  the  methods  used  to  analyze  and  visualize  model 
output  are  presented. 

3.2  IFM  Density  Backgrounds 

The  first  objective  was  to  obtain  IFM  electron  density  background  maps.  These 
served  two  purposes:  1)  provide  the  initial,  background  ionosphere  specification  for 
the  GAIM-GM  model;  and  2)  provide  a  background  observation  onto  which  a  syn¬ 
thetic  SED  could  be  imposed.  The  day  chosen  for  the  IFM  backgrounds  was  21  June 
2001  (day  172  of  the  year),  which  corresponded  to  the  summer  solstice  during  the 
peak  of  solar  cycle  23.  This  was  a  logical  choice  because  geomagnetic  storming  occurs 
most  often  during  solar  maximum.  Slant  TEC  measurements  of  the  IFM  backgrounds 
with  the  imposed  SED  could  then  be  cut  between  GPS  satellites  and  several  different 
distribution  of  GPS  ground  stations.  These  slant  TEC  measurements  were  the  only 
data  ingested  by  GAIM-GM. 

When  running  GAIM-GM,  two  days  of  IFM  output  are  needed  to  produce  a  one- 
day  GAIM-GM  specification  if  output  from  a  previous  IFM  run  is  not  available.  This 
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is  because  the  IFM  is  initiated  using  the  empirical  International  Reference  Ionosphere 
model  (IRI)  ( Space  Environment  Corporation ,  2002).  In  such  cases  the  IFM  output 
is  degraded,  so  it  is  necessary  to  run  IFM  for  two  warm-up  days,  followed  by  two 
good  days  in  order  to  arrive  at  a  steady-state  specification  of  the  ionosphere.  The 
IFM  input  parameters  for  the  two  good  days  (171  and  172)  are  listed  in  Table  2.  The 
input  parameters  for  the  two  warm-up  days  (169  and  170)  used  the  same  parameters 
as  day  171.  Note  that  only  the  Kp  and  Ap  indices  changed  between  days  171  and  172. 
The  F10.7  (10.7  cm  wavelength  solar  radio  flux)  and  F  10.1  A  (90-day  average  F10.7) 
indices  were  taken  from  historical  observations  and  reflect  solar  maximum  conditions 
during  early  2001. 


Table  2.  IFM  input  parameters 


Parameter 

Day  before  storm 

Day  of  storm 

Year 

2001 

2001 

Day 

171 

172 

Kp 

3 

6 

Ap 

14 

83 

FI  0.7 

205 

205 

FI  0.7 A 

153 

153 

A  Kp  of  3  was  used  on  day  171  to  reflect  a  relatively  quiet,  pre-storm  ionosphere 
state.  As  the  Kp  index  increased  to  6  on  day  172,  the  IFM  output  gave  significantly 
higher  TEC  values,  as  shown  in  Figure  10.  While  higher  TEC  is  often  seen  during 
geomagnetic  storming,  IFM  TEC  associated  with  elevated  Kp  indices  tends  to  be 
erroneously  high  due  to  missing  physics  in  the  IFM.  With  enough  assimilated  data, 
GAIM-GM  is  able  to  overcome  this  shortcoming  in  the  IFM  background  ( Schunk 
et  al. ,  2010). 
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Figure  10.  IFM  background  TEC  maps  at  1700  UT  for  days  171  (a)  and  172  (b)  of 
year  2001.  Much  larger  TEC  values  were  modeled  by  IFM  on  day  172  due  to  higher 
Kp  levels.  Unlike  Figure  7,  contour  smoothing  was  applied  to  these  images  to  remove 
the  pixelated  appearance 


3.3  Modeling  the  SED 

Once  the  IFM  runs  for  days  171  and  172  were  completed,  the  output  was  down¬ 
loaded  from  USU  servers  via  Secure  File  Transfer  Protocol  (SFTP)  as  Network  Com¬ 
mon  Data  Form  (NetCDF)  hies.  These  hies  were  then  opened  and  manipulated  in 
Matlab®.  An  IFM  NetCDF  hie  contains  one  simulated  day’s  worth  of  output  in  28 
different  variables.  Only  two  variables  were  used  during  the  process  of  modeling  the 
SED:  1)  three-dimensional  electron  density;  and  2)  TEC. 

The  entire  process  of  imposing  a  synthetic  SED  was  accomplished  in  the  Mat- 
LAB®utility  called  Storm,  m.  This  and  all  other  Matlab® program  used  during  this 
project  were  designed  by  the  author  at  the  Air  Force  Institute  of  Technology  (AFIT). 
Copies  of  these  programs  and  their  operating  instructions  may  be  requested  from  the 
author. 

Storm.m  hrst  made  two  adjustments  to  the  IFM  output  electron  density  array. 
The  hrst  adjustment  was  to  shave  exactly  3.75°  longitude  worth  of  data  from  both 
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edges  of  the  world  grid.  IFM  has  a  native  resolution  of  3°  latitude  by  7.5°  longitude 
and  gives  its  output  in  arrays  of  size  30  x  49.  The  1st  and  49th  longitude  indices 
(corresponding  to  0°  and  360°  longitude)  are  duplicates,  which  means  that  there  is 
an  extra  3.75°  of  longitude  at  the  eastern  and  western  edges  of  the  world  map.  These 
two  strips  of  data  ran  from  pole  to  pole  and  had  to  be  removed  so  that  the  IFM  would 
fit  exactly  within  the  boundaries  of  a  world  map. 

The  second  adjustment  was  to  increase  IFM’s  grid  resolution  to  match  GAIM- 
GM  in  regional  mode,  so  that  IFM  output  could  be  directly  compared  to  GAIM-GM 
output  at  a  later  point.  This  was  done  via  the  Matlab ®imresize  function  which  uses 
Lanczos  filtering,  a  simple  and  attractive  method  of  interpolating  between  sampled 
data  ( Duchon ,  1979).  The  Lanczos  filter  is  the  product  of  two  sine  functions,  where 
the  resulting  function  is  used  as  a  convolution  kernel  to  resample  the  input  field. 

An  example  of  IFM  TEC  output  before  and  after  interpolation  is  shown  in  Fig¬ 
ure  11.  The  IFM  output  in  11(a)  is  in  the  native  resolution  while  the  output  in  11(b) 
has  been  interpolated  and  resized  to  match  GAIM-GM’s  resolution.  Table  3  compares 
the  scale  and  resolution  of  IFM  global  output  to  GAIM-GM  regional  mode  output. 


Table  3.  IFM  output  size  and  resolution,  compared  to  GAIM-GM  regional  mode  size 
and  resolution.  The  IFM  had  to  be  resized  and  interpolated  to  match  GAIM-GM 
dimensions  so  that  direct  comparisons  over  CONUS  could  be  performed 


Parameter 

IFM 

GAIM-GM 

Latitude  Resolution 

3° 

1° 

Longitude  Resolution 

7.5° 

3.75° 

Western  Boundary 

0°  E 

225°  E 

Eastern  Boundary 

360°  E 

315°  E 

Northern  Boundary 

90°  N 

52°  N 

Southern  Boundary 

90°  S 

22°  N 

The  most  important  component  of  Storm. m  is  the  matrix  which  contains  the 
synthetic  SED  profile.  This  profile  was  patterned  after  the  behavior  of  real  SEDs, 
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TEC  Units  TEC  Units 


(a)  Native  resolution  of  3°  by  7.5° 


(b)  New  resolution  of  1°  by  3.75° 


Figure  11.  Interpolation  of  IFM  TEC  output  to  the  finer  GAIM-GM  resolution  for  year 
2001,  day  172,  time  1700  UT,  via  Lanczos  filtering  in  Matlab®.  Contour  smoothing 
was  not  applied  to  these  images,  hence  the  digitized  appearance 


several  of  which  have  been  well- documented  in  the  past  decade.  Ultimately  the  20 
November  2003  SED  analyzed  by  Foster  (2005)  was  chosen  as  the  basic  template 
for  the  synthetic  SED.  Figure  12  shows  GPS  TEC  observations  over  CONUS  at  the 
1945  UT  height  of  the  storm,  taken  from  about  450  North  American  GPS  sites. 
Additional  insight  on  the  SED’s  movement  was  gained  by  running  GAIM-GM  and 
ingesting  historical  slant  TEC  data  from  the  day  of  the  storm  to  replicate  the  SED 
features.  Using  89  North  American  GPS  sites,  GAIM-GM  produced  the  TEC  map 
and  electron  density  distribution  shown  in  Figure  13. 

For  simplicity  and  efficiency,  the  basic  SED  properties  were  captured  using  Fig¬ 
ure  12,  Figure  13,  and  Table  1.  Focus  was  placed  on  depicting  the  northwestward 
movement  of  a  plasma  parcel  within  the  broad  tongue  of  ionization,  which  in  Figure  12 
extends  from  the  U.S.  East  Coast  to  the  central  Canada.  The  finalized  synthetic  SED 
is  shown  in  Figure  14  in  terms  of  TEC  and  vertical  electron  density  profile. 

A  schematic  diagram  of  the  synthetic  SED’s  trajectory  is  given  in  Figure  15, 
identifying  its  precise  location  during  its  three-hour  transit  of  CONUS.  It  must  be 
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Figure  12.  GPS  vertical  TEC  observations  over  the  U.S.  at  1945  UT  during  the  height 
of  the  20  November  2003  geomagnetic  storm.  The  SED  is  the  tongue  of  ionization 
extending  from  the  Mid- Atlantic  coast  to  the  Great  Lakes.  450  GPS  receivers  provided 
vertical  TEC  data  to  create  this  image  (Adapted  from  Foster  (2005)) 
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Figure  13.  GAIM-GM  recreation  of  the  20  November  2003  geomagnetic  storm  at  1945 
UT.  Historical  slant  TEC  data  from  89  North  American  GPS  sites  (overlaid  as  black 
circles)  were  ingested  into  the  model,  (a)  shows  TEC  during  the  height  of  the  storm 
while  (b)  is  a  vertical  cross-section  of  electron  density,  taken  at  40°  N  latitude.  The 
model  resolution  at  the  time  was  1.333°  latitude  by  5°  longitude 
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Figure  14.  Synthetic  SED  in  an  IFM  background  at  1945  UT.  The  IFM  was  run  for 
year  2001,  day  172,  with  a  Kp  of  6  and  the  the  SED  was  superimposed  by  multiplying 
background  densities  between  200  —  600  km  by  a  factor  of  2.  (a)  shows  TEC  while  (b) 
is  a  vertical  cross-section  of  electron  density,  taken  at  40°  N  latitude.  The  model  grid 
resolution  is  1°  latitude  by  3.75°  longitude 


stressed  that  real  SEDs  consist  of  an  entire  plume  of  high  density  plasma  extending 
across  CONUS,  in  which  plasma  flows  northwestward  within  the  plume.  The  synthetic 
SED  in  this  research  is  a  much  more  basic  construct  that  only  models  a  parcel  of 
high  density  plasma  moving  along  the  path  that  the  plume  would  normally  outline. 
However,  the  properties  of  the  synthetic  SED  itself  closely  match  those  of  real  SEDs, 
as  can  be  verified  in  Table  4,  which  compares  the  properties  of  both  SED  types. 

Storm,  m  imposed  the  synthetic  SED  by  first  extracting  the  IFM  three-dimensional 
electron  density  output  at  each  time  step  for  the  target  day  (year  2001,  day  172) 
from  the  NetCDF  hie  and  loading  it  into  a  four- dimensional  read-only  array.  This 
array  consisted  of  latitude,  longitude,  altitude,  and  time  dimensions.  Storm.m  then 
multiplied  all  background  electron  densities  found  within  the  volumetric  boundary 
of  the  SED  region  by  a  factor  of  2  and  wrote  them  to  a  new  four-dimensional  storm 
density  array.  At  the  next  time  step,  the  SED  boundaries  were  propagated  exactly  2° 
latitude  to  the  north  and  3.75°  longitude  to  the  west  and  the  IFM  background  electron 
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Figure  15.  Trajectory  of  synthetic  SED  across  CONUS.  The  feature  originates  at  22°  N 
over  the  Atlantic  Ocean  at  1800  UT  and  travels  in  a  northwestward  direction  at  about 
0.5  km/s  until  disappearing  at  2100  UT  at  52°  N  over  lower  Manitoba,  Canada 


Table  4.  Synthetic  vs.  real  SED  properties 


Property 


Synthetic  SED 


Real  SED 


3  hours 


2-5  hours 


Duration 
Peak  time 
HmF2 

Base  elevation 

Ceiling  elevation 

Altitude  extent 

Length 

Width 

Velocity 

NmF2 

TEC  increase 


1800  -  2100  UT 

400  km 

200  km 

600  km 

400  km 

1500  km 

700  km 

0.5  km/s 

1  -  4  x  1012  m~3 

1.7  times 


2 

1730  -  2000  UT 

>  500  km 
Variable 
Variable 

300  -  500  km 

>  3000  km 
600  -  1000  km 
0.5  -  1  km/s 

1  -  3  x  1012  m~3 
2-4  times 
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densities  within  the  SED  volume  (from  the  read-only  array)  were  again  doubled.  By 
taking  densities  from  the  read-only  array  and  writing  them  to  the  new  storm  density 
array,  densities  within  the  SED  were  prevented  from  being  doubled  repeatedly.  After 
the  SED  was  fully  inserted,  Storm. m  used  an  algorithm  developed  by  Capt  Janelle 
Jenniges  to  transform  the  73-layer,  three-dimensional  electron  density  arrays  into  new 
TEC  maps  for  each  time  step  ( Jenniges ,  2011).  The  TEC  maps  were  retained  for 
later  comparison  to  GAIM-GM  output  while  the  four-dimensional  electron  density 
arrays  served  as  backgrounds  through  which  slant  TEC  measurements  were  taken. 

3.4  Creating  GPS  Ground  Station  Distributions 

Having  obtained  the  IFM  specifications  and  imposed  the  synthetic  SED,  the  next 
step  was  to  acquire  several  distributions  of  GPS  ground  stations.  Two  types  of  ground 
station  grids  were  created:  1)  grids  of  real  GPS  ground  stations  taken  from  the 
Continuously  Operating  Reference  Station  (CORS)  network,  hereafter  referred  to 
as  CORS  grids;  and  2)  grids  of  evenly  spaced  and  ideally  located,  artificial  ground 
stations,  hereafter  referred  to  as  Ideal  grids. 

In  the  CORS  grids  there  were  five  different  GPS  ground  station  distributions.  The 
first  distribution  was  taken  from  a  list  of  11  GPS  ground  stations  that  the  operational 
GAIM-GM  at  AFWA  had  routinely  ingested  during  February  2011.  It  was  found  that 
the  synthetic  SED’s  trajectory  put  the  AFWA  station  distribution  at  an  unreasonable 
disadvantage  for  detecting  and  replicating  the  storm  feature,  so  a  GPS  ground  station 
located  in  the  state  of  Georgia  was  added  to  the  AFWA  list,  bringing  the  total  number 
of  stations  to  12  (see  Table  5).  Justification  for  the  addition  of  this  12th  station  and 
its  impact  on  the  results  of  this  research  are  given  in  Appendix  B. 

Using  this  modified  12-station  AFWA  grid  as  the  baseline  distribution,  additional 
stations  were  added  in  increments  of  30,  100,  200,  and  400,  to  create  a  total  of  five 
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Table  5.  AFWA  GPS  stations 


Site 

ID 

Latitude  (°  N) 

Longitude  (°  E) 

Ellip.  height  (m) 

McDonald,  TX 

MDOl 

30.674 

255.991 

2005 

Brunswick,  GAa 

GABK 

31.150 

278.508 

-13 

Bermuda 

BRMU 

32.368 

295.308 

-10 

Pasadena,  CA 

JPLM 

34.202 

241.829 

424 

Pie  Town,  NM 

PIE1 

34.300 

251.881 

2348 

Goldstone,  CA 

GOLD 

35.421 

243.113 

987 

Washington,  DC 

USN3 

38.920 

282.930 

58 

Greenbelt,  MD 

GODE 

39.019 

283.177 

15 

North  Liberty,  I A 

NLIB 

41.769 

268.428 

208 

St.  John,  Canada 

STJO 

47.600 

307.320 

153 

Brewster,  WA 

BREW 

48.125 

240.323 

239 

Holberg,  Canada 

HOLB 

50.640 

231.870 

560 

a  Brunswick  was  not  part  of  the  original  11-station  AFWA  list.  Appendix  B 
explains  why  it  was  added  as  the  12th  station 


CORS  grids.  The  first  30  stations  added  to  the  AFWA  grid  were  hand-picked  to 
maximize  even  placement,  while  a  MATLAB®utility  called  StationS ort.m  selected  the 
100,  200,  and  400  additional  stations  via  a  special  filter  that  maximized  uniform 
spacing.  Figure  16  shows  the  locations  of  GPS  ground  stations  for  each  CORS  grid. 
It  is  important  to  emphasize  that  the  each  grid  builds  upon  the  previous  one,  such 
that  only  412  stations  were  used  overall. 

Attaining  a  perfectly  uniform  distribution  was  impossible  for  any  CORS  grid  be¬ 
cause  of  limitations  in  the  CORS  network  layout.  All  ^1900  CORS  occupy  only  259 
of  the  GAIM-GM  regional  mode  grid  squares  (there  are  720  total  grid  squares  in  the 
model),  creating  large  data  gaps  over  oceans  and  highly  saturated  areas  of  ground 
stations  over  land.  Although  the  AFWA+400  grid  was  distributed  as  evenly  as  possi¬ 
ble,  59%  of  station-hosting  model  grid  squares  contain  two  GPS  ground  stations  and 
one  square  contains  three  (see  Figure  17).  This  is  an  important  point  since  closely 
spaced  GPS  receivers  are  problematic  in  GAIM-GM,  as  TEC  measurements  become 
highly-correlated  and  are  overly  trusted  by  the  Kalman  filter  ( Thompson  et  al. ,  2005). 
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Latitude  Latitude 


225  240  255  270  285  300  315 

Longitude 


Longitude 

(b)  AFWA+30 


Longitude 

(c)  AFWA+100 


Longitude  Longitude 

(d)  AFWA+200  (e)  AFWA+400 

Figure  16.  CORs  grids  of  GPS  ground  stations.  These  grids  use  actual  station  locations 
from  the  CORS  network 
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225  240  255  270  285  300  315 

Longitude 


0  12  3 

Stations  Per  Grid  Square 


Figure  17.  GPS  ground  station  density  per  model  grid  square  for  the  AFWA+400 
grid.  106  model  grid  squares  have  only  one  station,  152  model  grid  squares  have 
two  stations,  and  one  model  grid  square  over  southern  California  has  three  stations. 
GAIM-GM  resolution  was  set  to  1°  latitude  by  3.75°  longitude 


The  Ideal  grids  were  designed  to  test  how  geo- location  accuracy  over  CONUS  im¬ 
proves  when  optimal  GPS  ground  station  spacing  is  applied,  completely  unrestricted 
by  oceans  and  terrain.  Figure  18  shows  the  four  Ideal  grids,  which  were  created  using 
21,  90,  360,  and  720  artificial  ground  stations.  Each  grid  builds  upon  the  last  in  the 
same  manner  as  the  CORS  grids.  The  main  differences  between  the  Ideal  grids  and 
the  CORS  grids  are:  1)  the  Ideal  grids’  ground  stations  are  artificial;  and  2)  the  Ideal 
grids  have  only  one  station  per  model  grid  square. 

For  both  the  CORS  and  Ideal  grids,  the  average  and  median  distances  separating 
nearest  neighbor  stations  are  listed  in  Table  6.  The  median  distance  is  a  more  rea¬ 
sonable  value  because,  unlike  the  average  distance,  it  is  not  skewed  by  contributions 
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Longitude 


(a)  Ideal  21 


Longitude 

(c)  Ideal  360 


Longitude 


(b)  Ideal  90 


225  240  255  270  285  300  315 

Longitude 


(d)  Ideal  720 


Figure  18.  Ideal  grids  of  GPS  ground  stations.  These  grids  represent  optimized  distri¬ 
butions  of  ground  stations  and  do  not  exist  in  reality 


from  geographically  isolated  stations  such  as  Bermuda.  Note  that  even  though  Ideal 
720  grid  has  more  stations  than  the  AFWA+400  grid,  the  distances  betweens  stations 
are  still  greater  because  the  stations  are  able  to  spread  out  over  the  oceans. 

3.5  Taking  Slant  TEC  Measurements 

Slant  TEC  measurements  are  derived  from  phase  differences  introduced  into  a  GPS 
satellite’s  signal  as  it  penetrates  the  ionosphere.  Unlike  simple,  single- channel  GPS 
devices  sold  commercially,  ground  receivers  at  official  GPS  monitoring  stations  are 
dual-channel  and  can  receive  both  the  LI  and  L2  signals  broadcast  by  GPS  satellites. 
Precise  measurement  of  the  time  delay  difference  between  these  two  signals  allows 
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Table  6.  Distance  between  GPS  ground  station  nearest  neighbors 


Distribution 

#  Stations 

Average  (km) 

Median  (km) 

Standard  Dev.  (km) 

AFWA 

12 

1580 

1321 

720 

AFWA+30 

42 

749 

686 

275 

AFWA+100 

112 

418 

362 

213 

AFWA+200 

212 

289 

252 

166 

AFWA+400 

412 

198 

171 

134 

Ideal  21 

21 

1325 

1250 

204 

Ideal  90 

90 

659 

621 

135 

Ideal  360 

360 

324 

308 

90 

Ideal  720 

720 

257 

247 

86 

the  TEC  along  the  signal  slant  path  to  be  calculated  (see  Appendix  A  for  further 
details). 

The  GPS  satellites  used  in  this  research  were  simulated,  but  their  time-dependent 
position  came  from  actual  archived  GPS  ephemerides.  These  archived  data  sets  from 
the  past  decade  to  the  present  are  available  from  NASA’s  Crustal  Dynamics  Data 
Information  System  (CDDIS)  in  the  form  of  .sp3  hies  and  can  be  downloaded  via  FTP 
( Crustal  Dynamics  Data  Information  System ,  2011).  GPS  constellation  ephemerides 
from  2  October  2010  were  selected  to  represent  the  GPS  constellation,  making  this 
research  more  relative  to  present  day  capabilities.  A  useful  diagram  of  these  satellites’ 
positions  over  the  course  of  a  24-hour  period  is  presented  in  Appendix  A,  Figure  48. 

With  the  GPS  ground  stations,  IFM  background,  and  simulated  GPS  constellation 
in  place,  the  slant  TEC  measurements  were  taken  using  a  MATLAB®utility  called 
SlantTEC.m.  This  program  required  three  data  sets  in  order  to  operate:  1)  the 
time-dependent  IFM  density  background  with  the  imposed  SED;  2)  a  list  of  GPS 
ground  station  coordinates  and  elliptical  elevations;  and  3)  a  .sp3  hie  containing  GPS 
satellite  ephemerides.  SlantTEC.m  used  a  15-minute  sampling  cadence  in  order  to 
match  GAIM-GM’s  data  assimilation  rate,  as  well  as  a  15°  Elevation  Mask  so  that 
slant  TEC  measurement  were  only  taken  from  satellites  at  least  15°  above  the  horizon 
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in  a  GPS  station’s  field  of  view.  The  output  from  SlantTEC.m  was  one  specially 
formatted  text  file  per  station,  per  simulated  hour,  that  could  be  uploaded  directly 
to  USU  servers  and  assimilated  by  GAIM-GM.  To  put  this  into  perspective,  using 
720  GPS  ground  stations  meant  that  17,280  text  hies  were  produced  to  capture  slant 
TEC  measurements  for  one  day.  These  hies  were  then  uploaded  to  GAIM-GM  via 
SFTP. 

The  portion  of  SlantTEC.m  that  actually  calculated  the  slant  TEC  was  adapted 
from  a  program  written  by  USU.  A  copy  of  USU’s  original  slant  TEC  calculation  code 
in  FORTRAN  can  be  obtained  by  contacting  the  Center  for  Atmospheric  and  Space 
Sciences.  Their  code  derives  the  TEC  between  two  points  given  in  Earth-centered, 
Earth-hxed  coordinates  by  dividing  the  slant  path  into  hxed  increments  and  linearly 
interpolating  the  electron  densities.  Each  interpolated  density  is  multiplied  by  the 
hxed  incremental  distance  in  order  to  determine  the  TEC  for  the  slant  path  segment. 
The  TEC  values  are  then  totaled  to  produce  the  slant  TEC  value  for  the  entire  slant 
path  (Nava,  2011). 

To  ensure  the  accuracy  of  the  USU  slant  TEC  code  and  quantify  the  error  its 
density  interpolation  scheme  introduced,  a  simple  validation  was  performed.  In  Mat- 
LAB®,  720  simulated,  non-moving  GPS  satellites  were  held  hxed  over  an  interpolated 
IFM  density  background’s  720  model  grid  squares.  These  simulated  satellites  then 
took  vertical  TEC  measurements  at  each  15-minute  time  step  for  all  24  hours  of  IFM 
year  2001,  day  172  output  (without  the  imposed  SED).  The  vertical  slant  TEC  pierc¬ 
ings  from  every  grid  square  were  then  compared  to  actual  IFM  TEC  values,  as  taken 
directly  from  the  IFM’s  NetCDF  output  hies.  Figure  19  shows  the  average  difference 
for  each  grid  square  over  the  24-hour  period.  Dark  shades  of  blue  and  red  mark  lo¬ 
cations  where  the  largest  differences  (i.e.  errors)  occurred,  while  white  grid  squares 
indicate  locations  with  very  small  errors.  The  average  TEC  difference  for  the  entire 
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-3-2-10123 

Percent 

Figure  19.  Average  vertical  slant  TEC  error  per  grid  square  over  a  24-hour  period, 
taken  for  year  2001,  day  172.  Vertical  piercings  of  the  IFM  density  background  were 
performed  at  every  15-minute  time  step  and  compared  to  actual  TEC  values.  The 
mean  absolute  difference  for  all  vertical  slant  TEC  measurements  was  2.7% 


map  was  2.7%  while  the  maximum  difference  of  any  single  grid  square,  for  all  time 
steps,  was  3.2%.  The  greatest  slant  TEC  errors  generally  occurred  in  areas  of  strong 
TEC  gradients,  such  as  along  the  northern  edge  of  the  equatorial  anomaly. 

The  average  uncertainty  in  GAIM-GM  output  was  5.7%,  based  solely  on  the  3.0 
TEC  unit  Phase  TEC  error  assigned  to  each  slant  TEC  data.  Although  there  was  no 
actual  Phase  TEC  error  in  this  research  because  all  measurements  were  simulated, 
there  were  two  other  errors  types  that  needed  to  be  accounted  for  by  GAIM-GM  (see 
Table  7).  For  simplification  purposes,  the  Phase  TEC  error  was  retained  as  a  proxy 
for  the  other  two  error  types,  which  were  due  to  1)  Slant  TEC  pierce  point  averaging; 
and  2)  IFM/GAIM-GM  ceiling  difference. 

Slant  TEC  pierce  point  averaging  error  caused  an  average  of  2.7%  uncertainty  and 
occurred  as  SlantTEC.m  took  the  mean  density  of  all  model  grid  squares  adjacent  to 
the  pierce  grid  square,  at  every  path  increment.  This  localized  interpolation  helped 
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Table  7.  Average  uncertainty  in  GAIM-GM  output  due  to  slant  TEC  measurements. 
The  sub-types  stem  from  methods  used  in  this  research  and  are  accounted  for  by  the 
Phase  TEC  error 


Type 

Amount 

Explanation 

Phase  TEC  error 

5.7% 

Derived  from  the  3.0  TEC  unit  error  assigned  to 
every  synthetic  slant  TEC  measurement  (a  typical 
value  for  real  slant  TEC) 

Sub- type 

Amount 

Explanation 

Slant  TEC  pierce 
point  averaging 

2.7% 

See  Figure  19.  USU  slant  TEC  code  takes  the 
mean  density  of  all  grid  squares  adjacent  to  the 
pierced  grid  square,  for  every  path  increment.  This 
is  done  to  represent  the  continuous  distribution  of 
densities  in  the  real  ionosphere 

IFM/GAIM-GM 
ceiling  difference 

1.4% 

Average  amount  of  TEC  that  the  1400  -  1600  km 
layer  of  the  IFM  background  contributes  to  slant 
TEC  measurements.  GAIM-GM  only  has  a  ceiling 
of  1400  km  but  still  assimilates  the  full  slant  TEC 
value  taken  from  the  IFM  background 

represent  the  continuous  distribution  of  densities  in  the  real  ionosphere  -  unlike  the 
discontinuous,  pixelated  densities  in  the  IFM.  The  IFM/GAIM-GM  ceiling  differ¬ 
ence  error  produced  an  average  of  1.4%  uncertainty.  This  occurred  because  the  IFM 
extends  200  km  higher  than  GAIM-GM  and  extra  topside  ionosphere  plasma  was 
integrated  into  slant  TEC  measurements.  GAIM-GM  accepted  this  extra  fraction  of 
TEC  and  applied  it  towards  its  90  -  1400  km  model  region.  Combined,  these  two 
error  sources  did  not  amount  to  the  5.7%  reported  Phase  TEC  error,  but  the  higher 
5.7%  value  was  nonetheless  kept  because  it  represents  the  amount  of  error  expected 
in  real  world  slant  TEC  data. 

3.6  Running  GAIM-GM 

GAIM-GM  is  a  highly  sophisticated  computer  model  and  detailed  instructions  re¬ 
garding  its  operation  are  contained  in  the  User’s  Guide,  which  may  be  requested  from 
USU.  However,  the  basic  process  for  using  GAIM-GM  to  assimilate  slant  TEC  and 


37 


replicate  the  synthetic  SED  are  shown  in  a  flowchart  in  Figure  20.  The  shaded  right 
half  of  the  flowchart  depicts  the  operations  performed  physically  at  AFIT,  while  the 
unshaded  left  half  of  the  flowchart  depicts  operations  that  occurred  at  USU.  Arrows 
flowing  across  shaded  regions  represent  SFTP  actions.  Dark  grey  ovals  represent 
different  inputs  and  light  grey  boxes  depict  the  computer  programs  which  utilized 
those  inputs.  After  all  data  were  assimilated  and  GAIM-GM  produced  its  solution 
in  NetCDF  hie  format,  those  output  hies  were  downloaded  and  a  variety  of  Mat- 
LAB®analysis  utilities  processed  the  raw  numeric  results  into  useful  hgures  and  plots. 

A  total  of  10  cases  were  performed  under  varying  model  conhgurations,  SED 
characteristics,  and  GPS  ground  station  distributions,  for  26  runs  in  all.  Table  8 
gives  a  complete  list  of  these  runs  and  their  variations.  More  details  on  each  run  will 
be  presented  in  the  next  chapter. 

The  GAIM-GM  model  used  in  this  research  was  based  on  version  2.8.1  that  AFWA 
uses  operationally,  but  with  three  slight  adjustments.  First,  the  GPS  ground  station 
and  satellite  biases  were  deactivated  because  the  simulated  GPS  ground  stations  and 
satellite  constellation  had  no  biases.  Second,  GAIM-GM’s  internal  algorithm  that 
subtracts  out  the  plasmaspheric  contribution  to  TEC  was  deactivated.  Aside  from 
the  minor  GAIM-GM/IFM  ceiling  difference  described  in  Table  7,  the  IFM  density 
background  had  precisely  known  values  and  there  was  no  extra  plasmasphere  layer. 
Third,  Cases  2  through  7  were  run  using  a  special  version  of  GAIM-GM  that  assigned  a 
high  amount  of  trust,  or  heavier  weight,  to  assimilated  data.  While  this  was  useful  for 
testing  GAIM-GM’s  capabilities  to  resolve  a  relatively  small  storm  feature,  there  are 
hazards  in  placing  too  much  trust  in  observations,  especially  if  the  data  is  erroneous. 
Case  6  illustrates  this  point  by  reducing  all  slant  TEC  measurements  at  the  Freeport, 
Ohio  GPS  station  by  a  factor  of  4  to  simulate  bad  data,  which  in  reality  could  be 
caused  by  instrument  malfunction  or  improper  calibration. 
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Figure  20.  GAIM-GM  model  run  flowchart 
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AFIT 


3.7  Analysis  Techniques 


Several  analysis  methods  were  applied  to  extract  useful  results  from  the  GAIM- 
GM  output.  The  statistics  derived  from  these  methods  were  expressed  in:  1)  dif¬ 
ferences  (for  TEC  and  electron  density);  2)  percentage  differences;  3)  mean  absolute 
errors  (MAEs);  and  4)  skill  scores. 

Difference  plots  are  self-explanatory  and  provide  a  fast  way  to  visualize  the  dif¬ 
ferences  between  two  model  runs.  Two-dimensional  latitude  by  longitude  plots  were 
used  to  compare  TEC,  while  two-dimensional  longitude  by  altitude  electron  density 
plots  were  used  to  identify  differences  in  electron  density  vertical  profiles.  The  dif¬ 
ference  at  a  single  grid  square  is  simply  the  value  of  first  output  minus  the  second 
output.  The  percentage  difference  at  a  single  grid  square,  DIFFm ,  is: 


DIFFm 


( SPECIm  -  REFm\ 
V  SPECIm  J 


x  100 


(4) 


where  SPECIm  is  the  GAIM-GM  specification  of  interest  and  REFm  is  a  reference 
that  can  be  either  a  GAIM-GM  or  IFM  specification. 

MAE  is  an  accuracy  measure  that  is  sometimes  used  instead  of  root  mean  square 
error  (RMSE)  because  it  is  not  as  sensitive  to  outlying  data  ( Hyndman  and  Koehler , 
2005).  For  these  results,  MAE  applied  only  to  TEC  output  and  was  used  to  measure 
how  much  the  entire  GAIM-GM  map  differed  from  the  IFM  background  at  a  given 
time  step.  One  value  is  generated  for  each  time  step,  according  to  the  formula: 


MAE  =  jf  E  I TEC GAM  -  TECifm 


(5) 


where  TECifm  is  the  IFM  TEC  value  and  TECq^im  is  the  GAIM-GM  TEC  value 
for  a  given  GPS  ground  station  grid. 

Skill  scores  are  commonly  used  by  meteorologists  but  can  easily  be  applied  to  space 
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weather  applications.  These  scores  are  generally  defined  as  measures  of  the  relative 
accuracy  of  forecasts  produced  by  two  forecasting  systems,  one  of  which  is  a  reference 
system  -  such  as  climatology.  Positive  skill  (in  other  words,  favorable  differences  in 
accuracy)  typically  represents  a  minimal  level  of  acceptable  performance  for  a  set  of 
forecasts.  A  forecast  receives  a  positive  skill  score  if  it  matches  observations  more 
accurately  than  the  climatology;  if  it  is  less  accurate  than  the  climatology,  it  receives 
a  negative  skill  score  ( Murphy ,  1988). 

In  the  context  of  this  research,  skill  scores  measure  how  well  a  GAIM-GM  model 
specification  reproduces  the  IFM  background  with  imposed  SED,  compared  to  a  refer¬ 
ence  GAIM-GM  specification  that  is  also  attempting  to  reproduce  the  IFM.  The  ref¬ 
erence  specification  is  the  GAIM-GM  AFWA  grid  solution,  which  represents  AFWA’s 
current  capabilities  using  its  12  GPS  ground  stations  over  the  vicinity  of  CONUS. 
The  formula  for  these  skill  scores  was  taken  from  a  previous  study  that  evaluated 
how  well  GAIM-GM  and  the  IFM  performed  with  respect  to  the  IRI  over  Australia 
during  2004  ( Sojka  et  al. ,  2007). 

To  find  the  skill  score,  first  the  RMSE  scores  are  determined  for  the  desired  GAIM- 
GM  grid  specification  and  also  for  the  GAIM-GM  AFWA  specification,  given  respec¬ 
tively  as: 

Score  GAIMGrid  =  Cc  (e 

and 

Score  GAIMafwa  =  y/T,  ^ecGm  -  tec™)2  (7) 

These  scores  are  then  combined  to  generate  the  GAIM  skill  score ,  which  is: 


GAIM$«fLL  = 


Score  GAIMGrid  \ 

1.0  -  — - -  x  100 

Score  GAIMafwa  J 


(8) 


The  skill  scores  are  expressed  as  points  and  are  positively  bounded  at  100  when 
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a  GAIM-GM  grid  specification  exactly  matches  the  IFM  background.  Positive  skill 
scores  imply  that  the  grid  specification  was  more  accurate  than  the  reference  AFWA 
specification,  while  negative  scores  indicate  the  AFWA  specification  was  more  accu¬ 
rate.  Values  of  zero  imply  that  both  specifications  were  equivalent  in  reproducing  the 
IFM  background. 

Although  skill  scores  and  MAEs  are  typically  taken  over  an  entire  map,  they 
can  also  be  computed  for  just  certain  portions  of  the  map.  With  careful  coordinate 
indexing  they  can  be  customized  to  track  a  moving  feature,  or  conversely  everything 
except  the  moving  feature,  which  is  exactly  what  the  MATLAB®utility  ScoreGAIM.m 
did.  It  took  the  SED’s  exact  size  and  trajectory  from  matrices  in  Storm. m  and  created 
a  sampling  grid  that  moved  with  the  SED  (see  Figure  15)  and  calculated  skill  scores 
and  MAEs  for  just  that  area.  It  also  derived  skill  scores  and  MAEs  for  everything 
except  the  SED  by  treating  the  moving  sampling  grid  as  a  void  in  the  map  and  scoring 
everything  but  that  void.  Developing  statistics  for  these  three  domains  (entire  grid, 
the  SED  only,  and  everything  but  the  SED)  allowed  for  a  more  robust  analysis  of 
GAIM-GM’s  accuracy. 

In  addition  to  the  skill  score  and  MAE  plots,  unique  skill  score  maps  were  created 
to  display  scores  over  a  two-dimensional  latitude  by  longitude  grid.  The  formula  for 
each  grid  square  on  these  maps  is  similar  to  Equation  8  except  averages  are  not  taken 
over  the  entire  grid.  Each  grid  square  is  scored  as: 


n  AT  A/tGrid 
Or/ii  1V1 SKILL  MAP 


/  (TECIFM-TEC%tfM)2\ 

V  (TECIFM-TEC£FAyMA?) 


x  100 


(9) 


As  shown  in  the  next  chapter,  these  skill  score  maps  use  a  64-hue  color  scale, 
where  green  tones  show  positive  scores  and  brown  tones  show  negative  scores.  White 
areas  represent  zero  values  where  the  grid  specification  was  equivalent  to  the  AFWA 
specification. 
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3.8  Data  Visualization 


Just  as  skill  score  maps  can  provide  more  insight  than  numbers  alone,  data  visu¬ 
alization  allows  for  quick  and  easy  interpretation  of  the  enormous  amounts  of  data 
produced  by  GAIM-GM  during  this  research.  All  data  visualization  products  were 
created  in  Matlab®,  to  include:  1)  GPS  ground  station  maps;  2)  GPS  ground  sta¬ 
tion  density  maps;  3)  TEC  maps;  4)  TEC  difference  maps;  5)  electron  density  latitude 
slices;  6)  electron  density  difference  profiles;  7)  MAE  plots;  8)  skill  score  graphs;  and 
9)  skill  score  maps. 

The  GPS  ground  station  and  density  maps  were  relatively  easy  to  create,  requiring 
nothing  but  a  spreadsheet  of  4-digit  station  identifiers  and  coordinates  as  input.  The 
TEC  and  electron  density  maps  were  a  more  difficult  chore  to  produce  and  required 
extracting  thousands  of  two-dimensional  and  three-dimensional  arrays  from  NetCDF 
output  files  and  loading  them  into  four- dimensional  arrays  under  a  stringent  system 
of  bookkeeping.  Although  massive  in  size,  the  four- dimensional  arrays  dramatically 
reduced  the  number  of  variables  required  to  contain  the  24  hours  worth  of  TEC  and 
electron  density  specifications  for  each  individual  model  run.  The  tradeoff  for  ease 
of  use  was  higher  computational  demands.  Workstations  lacking  64-bit  operating 
systems  and  at  least  4  gigabytes  of  random  access  memory  were  incapable  of  running 
some  of  the  MATLAB®utilities.  Copies  of  the  MATLAB®utilities  from  this  project 
with  concise  operating  instructions  instructions  may  be  requested  from  the  author. 
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Table  8.  Complete  list  of  GAIM-GM  regional  mode  runs 


Case 

GAIM-GM 

GPS  Grid 

IFM 

Variations 

0a 

Normal 

None 

Normal 

1-A 

Normal 

AFWA 

Normal 

1-B 

Normal 

AFWA+30 

Normal 

1-C 

Normal 

AFWA+100 

Normal 

1-D 

Normal 

AFWA+200 

Normal 

1-E 

Normal 

AFWA+400 

Normal 

2-Ab 

Modified 

AFWA 

Normal 

2-B 

Modified 

AFWA+30 

Normal 

2-C 

Modified 

AFWA+100 

Normal 

2-D 

Modified 

AFWA+200 

Normal 

2-E 

Modified 

AFWA+400 

Normal 

3-A 

Modified 

Ideal  21 

Normal 

3-B 

Modified 

Ideal  90 

Normal 

3-C 

Modified 

Ideal  360 

Normal 

3-D 

Modified 

Ideal  720 

Normal 

4 

Modified 

AFWA+100 

Normal 

SED  raised  up  in  altitude 

5 

Modified 

AFWA+400 

Normal 

Time  Constant  reduced  to  1 

6 

Modified 

AFWA+100 

Normal 

One  GPS  station  gave  bad  data 

7 

Modified 

Ideal  21 

Normal 

GPS  station  grid  reversed  left-right 

8-Ac 

Normal 

AFWA 

Depleted 

8-B 

Normal 

AFWA+30 

Depleted 

8-C 

Normal 

AFWA+100 

Depleted 

8-D 

Normal 

AFWA+200 

Depleted 

8-E 

Normal 

AFWA+400 

Depleted 

9 

Normal 

AFWA 

Depleted 

Time  Constant  reduced  to  1 

10 

Normal 

AFWA 

Depleted 

Elevation  Mask  raised  to  45° 

a  Zero  station  baseline  run.  Since  there  was  no  data  assimilation,  the  version  of 
GAIM-GM  used  and  GPS  station  grid  type  had  no  effect  on  the  output 
b  Cases  2- A  through  7  were  not  run  in  GAIM-GM  2.8.1,  but  a  modified  version  of 
GAIM-GM  that  assigns  more  trust  to  assimilated  data 

c  Cases  8- A  through  10  took  all  slant  TEC  measurements  though  an  IFM  background 
where  all  densities  north  of  30°  N  latitude  (excluding  the  SED)  were  depleted  by  a 
factor  of  4  during  the  storm  period,  1800  -  2045  UT 
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IV.  Results  and  Analysis 


4.1  Chapter  Overview 

The  objective  of  this  research  is  to  quantify  how  many  ground  GPS  receivers 
are  needed  in  CONUS  to  provide  GAIM-GM  with  sufficient  data  to  create  reliable 
GPS/TEC  correction  maps  during  geomagnetic  storms.  The  results  presented  in  this 
chapter  will  show  how  varying  the  number  of  GPS  ground  stations,  under  several 
different  scenarios,  affected  the  accuracy  of  TEC  maps.  Ten  separate  cases  were 
examined  and  each  case  comprises  one  section  of  this  chapter. 

Sections  1  and  2  discuss  using  the  CORS  grids  to  reproduce  the  synthetic  SED  in 
the  normal  and  modified  versions  of  GAIM-GM,  respectively,  while  Section  3  discusses 
applying  the  Ideal  grids  to  modified  GAIM-GM.  Sections  4  through  7  examine  the 
following  unique  situations:  4)  modeling  a  high- altitude  synthetic  SED;  5)  reducing 
the  Time  Constant  from  5  down  to  1;  6)  introducing  simulated  erroneous  data  into 
GAIM-GM  via  one  bad  station;  and  7)  reversing  the  layout  of  the  21-station  Ideal 
grid.  The  8th  section  discusses  the  results  of  depleting  the  IFM  background  around 
the  SED  to  simulate  the  real  SED  from  20  Nov  2003  (see  Figure  12).  Lastly,  Sections 
9  and  10  examine  the  effects  of  reducing  the  Time  Constant  from  5  to  1  and  raising 
the  Elevation  Mask  from  15°  to  45°  in  the  depleted  IFM  environment. 

4.2  Baseline  Run  with  Zero  GPS  Ground  Stations 

The  baseline  run  was  simply  the  default  output  from  GAIM-GM  when  no  slant 
TEC  measurements  were  assimilated.  It  is  important  to  understand  how  GAIM-GM 
performs  in  the  absence  of  data,  using  only  the  internal  IFM  background  and  error 
covariance  matrix  to  build  a  solution  (see  Section  2.5.3).  First  note  that  because  of 
the  error  covariance  matrix,  which  contains  an  IFM  climatology  from  1107  model 
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runs,  GAIM-GM  does  not  make  an  exact  replica  of  the  IFM  density  background  in 
the  absence  of  external  data.  Rather,  the  Kalman  filter  imposes  perturbations  from 
the  error  covariance  matrix  onto  the  IFM  background  densities,  yielding  a  slightly 
modified  solution. 

Figure  21  illustrates  this  by  comparing  IFM  TEC  output  to  GAIM-GM  TEC 
output  with  no  data  assimilation  for  year  2001,  day  172,  at  1945  UT.  The  IFM’s 
geophysical  input  parameters  were  listed  previously  in  Table  2  and  are  the  same 
for  every  model  run  in  this  project.  Note  in  21(b)  how  GAIM-GM  has  lower  TEC 
values  than  the  IFM,  by  an  average  of  3.8  TEC  units.  The  IFM  tends  to  produce 
erroneously  high  TEC  values  during  geomagnetic  activity  ( Schunk  et  al .,  2010),  and 
so  even  without  ingesting  data  GAIM-GM  provided  a  specification  with  slightly  lower 
TEC  values  than  the  IFM. 

4.3  Case  1:  Normal  GAIM-GM  and  CORS  Grids 

This  case  examined  GAIM-GM’s  ability  to  reproduce  the  synthetic  SED  using 
the  five  different  CORS  grids.  Figure  22  shows  the  IFM  TEC  background  compared 
to  the  GAIM-GM  TEC  specifications  at  1945  UT.  The  GAIM-GM  specifications 
visually  showed  increasing  improvement  in  matching  the  IFM  background  as  GPS 
ground  stations  were  added  to  the  existing  AFWA  grid,  due  to  more  slant  TEC 
measurements. 

The  vertical  electron  density  profiles  at  40°  N  latitude  for  each  model  run  also 
improved  as  the  station  count  increased.  Adding  30  ground  stations  to  the  AFWA  list 
caused  the  most  visible  improvement  from  one  run  to  the  next,  as  seen  in  Figure  23. 
The  unit  for  these  density  plots  is  electrons  /  cm~3 .  The  highest  density  obtained 
in  the  AFWA  grid  (see  Figure  23(a))  was  1.35xl06  electrons / cm~3 ,  compared  to 
1.69xl06  electrons / cm~3  in  the  AFWA+30  grid  (see  Figure  23(b)).  Clearly  more 
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Figure  21.  Comparison  of  IFM  and  GAIM-GM  TEC  specifications  using  zero  GPS 
ground  stations  at  1945  UT.  The  difference  of  (b)  -  (a)  is  given  in  (c).  Even  though 
GAIM-GM  falls  back  on  the  IFM  density  background  when  no  data  is  available,  GAIM- 
GM’s  internal  error  covariance  matrix  adjusts  the  solution  based  on  climatology  from 
1107  archived  IFM  runs  (see  Section  2.5.3),  causing  the  output  to  vary  from  the  IFM 
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Figure  22.  GAIM-GM  TEC  reproduction  of  the  synthetic  SED  at  1945  UT,  using  slant 
TEC  data  from  the  CORS  grids.  The  synthetic  SED  in  the  IFM  background  is  shown 
in  (a)  as  a  reference  and  the  locations  of  the  12  AFWA  ground  stations  are  shown  as 
black  double  circles  in  (b).  GAIM-GM  resolution  for  this  and  all  subsequent  figures  is 
1°  latitude  by  3.75°  longitude 
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Figure  23.  GAIM-GM  electron  density  reproduction  of  the  synthetic  SED,  using  the 
CORS  grid  at  1945  UT.  Latitude  slices  at  40°  N  are  shown  for  the  AFWA  grid  (a) 
and  AFWA+30  grid  (b).  The  high-density  core  of  the  SED  is  resolved  much  more 
accurately  when  30  stations  are  added  to  the  AFWA  grid 


slant  TEC  piercings  from  a  denser  ground  station  network  allowed  the  AFWA+30 
grid  to  take  more  cuts  through  the  SED  core,  improving  its  solution.  The  actual  IFM 
background  density  reached  as  high  as  2.29xl06  electrons / cm~3  in  the  core  of  the 
IFM  synthetic  SED  (not  shown).  Output  from  the  AFWA+100,  AFWA+200,  and 
AFWA+400  runs  appeared  very  similar  to  Figure  23(b)  and  are  not  presented. 

A  direct  comparison  of  the  GAIM-GM  TEC  to  the  IFM  TEC  for  the  AFWA 
and  AFWA+30  grids  is  shown  in  Figure  24.  This  type  of  plot  is  useful  for  quickly 
identifying  where  GAIM-GM  performed  best  at  reproducing  the  IFM  background. 
The  AFWA+30  grid  results  in  Figure  24(b)  have  much  more  white  areas,  indicating 
overall  better  IFM  reproduction  than  the  AFWA  grid  results.  Note  the  red  pixels 
(positive  differences)  immediately  outside  the  SED  area  in  Figure  24(b).  Here  GAIM- 
GM  had  difficulty  resolving  the  sharp  edge  of  the  SED  where  IFM  densities  abruptly 
doubled,  and  instead  smoothed  the  high  TEC  values  over  a  broader  area. 

The  percentage  TEC  difference  between  the  AFWA  and  AFWA+30  GAIM-GM 
runs  is  shown  in  Figure  25(a).  To  the  right,  Figure  25(b)  is  a  skill  score  map  showing 
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Figure  24.  TEC  difference  between  GAIM-GM  and  IFM  background  at  1945  UT.  The 
IFM  was  subtracted  from  the  GAIM-GM  AFWA  grid  specification  to  create  (a),  while 
(b)  shows  the  IFM  subtracted  from  the  GAIM-GM  AFWA+30  grid  specification.  The 
SED’s  location  is  clearly  indicated  by  the  dark  blue  patches  where  GAIM-GM  differed 
most  from  the  IFM.  White  regions  indicate  where  GAIM-GM  reproduced  the  IFM 
background  with  very  little  error 


where  the  AFWA+30  grid  improved  over  the  AFWA  grid  in  reproducing  the  IFM 
background.  It  is  important  to  understand  that  TEC  difference  plots  do  not  show 
accuracy  improvement,  whereas  a  skill  score  map  does.  Notice  how  the  AFWA+30 
grid  had  positive  skill  scores  directly  over  the  SED  but  negative  skill  scores  around  the 
SED  periphery,  due  to  the  sharp  gradient  effect  explained  in  the  previous  paragraph. 
Also  notice  how  adding  30  stations  to  CONUS  understandably  did  little  to  improve 
accuracy  over  the  data  sparse  Atlantic  Ocean,  as  indicated  by  the  white  area  in  the 
right  side  of  Figure  25(b). 

Average  skill  scores  per  time  step  for  all  five  runs  during  the  1800  -  2045  UT 
storm  period  are  presented  in  Figure  26.  A  plot  of  the  SED’s  trajectory  is  given  in 
Figure  26(a)  as  a  reference.  For  the  SED-only  domain,  all  scores  were  nearly  identical 
until  about  1900  UT  because  up  to  that  time  Bermuda  was  the  only  data  source  for 
all  grids.  As  the  SED  neared  the  East  Coast,  the  higher-density  grids  provided  more 
slant  TEC  data,  causing  their  skill  scores  to  increase  faster  than  the  low  density  grids. 
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Figure  25.  TEC  difference  and  skill  score  map  for  the  AFWA+30  grid  at  1945  UT,  com¬ 
pared  to  the  AFWA  grid.  The  AFWA  grid  TEC  was  subtracted  from  the  AFWA+30 
grid  TEC  to  yield  (a),  while  (b)  is  the  AFWA+30  grid  skill  score.  Greens  in  the  skill 
score  map  show  where  the  AFWA+30  grid  was  more  accurate  than  the  AFWA  grid  in 
reproducing  the  IFM  background,  while  brown  colors  show  where  it  was  less  accurate. 
Note  that  TEC  difference  maps  cannot  indicate  areas  of  accuracy  improvement  over 
the  AFWA  grid,  whereas  skill  score  maps  can 


A  different  phenomenon  was  responsible  for  the  skill  score  pattern  seen  in  Fig¬ 
ure  26(d),  where  scores  for  all  areas  outside  the  SED  declined  once  the  SED  reached 
the  mainland.  Recall  that  the  densely  populated  GPS  station  grids  provided  more 
slant  TEC  piercings  through  the  SED  than  the  sparsely  populated  grids.  Because 
GA1M-GM  could  not  resolve  the  sharp  boundary  of  the  SED,  it  spread  high-density 
plasma  perturbations  from  within  the  SED  to  areas  outside  of  the  SED.  These  errors 
accumulated  with  time  as  GAlM-GM’s  default  Time  Constant  setting  of  5  hours  held 
the  high-density  perturbations  in  the  grid  and  propagated  them  downstream  (i.e.  to 
the  west),  causing  skill  scores  outside  the  SED  to  decrease.  By  2045  UT  all  skill 
scores  had  begun  to  improve,  most  likely  because  SED  perturbations  were  less  dense 
at  higher  latitudes  and  caused  smaller  TEC  errors. 

Table  9  gives  an  average  of  the  time-dependent  skill  scores  for  each  run  over 
the  1800  -  2045  UT  storm  period,  ft  must  be  noted  that  the  skill  scores  for  all 
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Figure  26.  Skill  scores  for  the  CORS  grids  during  the  storm  period  from  1800  to  2045 
UT.  Compared  are  the  accuracy  improvement  over  the  AFWA  grid  for  three  domains: 
(b)  the  entire  map,  (c)  SED  only,  and  (d)  everything  outside  the  SED.  The  trajectory 
of  the  SED  is  shown  in  (a)  as  a  reference.  Scores  above  zero  indicate  the  grid  was  more 
accurate  at  reproducing  the  IFM  background  than  the  AFWA  grid,  while  scores  below 
zero  indicate  worse  accuracy 
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results  depended  on  the  SED  size  and  trajectory,  as  well  as  the  specific  combination 
of  settings  used  in  GAIM-GM.  The  zero  station  scores  were  all  negative,  indicating 
worse  accuracy  than  that  produced  by  the  12-station  AFWA  grid.  Skill  scores  for  the 
entire  map  stayed  within  1  point  for  all  five  runs,  while  SED-only  skill  scores  showed 
steady  improvement  as  station  counts  increased.  Conversely,  scores  for  areas  outside 
the  SED,  while  better  in  all  cases  than  the  AFWA  grid  scores,  actually  decreased 
as  the  number  of  stations  increased.  The  price  for  better  accuracy  in  the  SED-only 
domain  was  a  reduction  in  accuracy  everywhere  else,  but  most  notably  in  the  region 
immediately  surrounding  the  SED,  as  was  shown  in  Figure  25(b).  Because  GAIM- 
GM  could  not  resolve  the  SED’s  sharp  edges  well,  adding  more  GPS  stations  caused 
more  erroneous  density  perturbations  to  accumulate  around  the  SED  periphery. 


Table  9.  Summary  of  Case  1  average  skill  scores.  Standard  deviations  are  given  in 
parenthesis  to  the  right  of  each  score 


Domain 

No  stations 

AFWA 

+30 

AFWA 

+100 

AFWA 

+200 

AFWA 

+400 

Entire  map 
SED  only 
Outside  SED 

-19.0  (11.0) 
-10.2  (7.5) 
-39.1  (24.0) 

13.1  (5.1) 
13.5  (10.5) 

21.1  (13.1) 

13.0  (4.0) 
17.5  (14.4) 
18.4  (14.8) 

13.9  (4.6) 
20.8  (16.6) 
14.4  (15.7) 

14.0  (4.9) 
24.1  (16.9) 
6.7  (16.0) 

4.4  Case  2:  Modified  GAIM-GM  and  CORS  Grids 

The  same  CORS  grids  from  Case  1  were  applied  in  Case  2,  which  used  the  modified 
version  of  GAIM-GM  that  assigns  more  weight  to  ingested  data.  Figure  27  shows  the 
TEC  results  from  the  five  model  runs  at  1945  UT,  compared  to  the  IFM  background 
with  imposed  SED.  In  contrast  to  Figure  22,  the  SED  was  resolved  more  clearly. 
Notice  there  was  very  little  change  in  TEC  beyond  the  AFWA+30  grid. 

Comparing  GAIM-GM  minus  IFM  difference  plots  at  1945  UT  to  those  from 
Case  1  (see  Figure  24),  the  modified  GAIM-GM  solution  did  better  at  reproducing 
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Figure  27.  GAIM-GM  TEC  reproduction  of  the  synthetic  SED  at  1945  UT,  using 
GPS  slant  TEC  data  from  the  CORS  grid.  Note  that  the  synthetic  SED  in  the  IFM 
background  is  shown  in  (a)  as  a  reference  and  that  the  12  AFWA  ground  stations  are 
plotted  in  (b)  as  black  circles 
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Figure  28.  TEC  difference  between  GAIM-GM  and  IFM  background  at  1945  UT.  The 
IFM  was  subtracted  from  the  GAIM-GM  AFWA  grid  specification  to  create  (a),  while 
(b)  shows  the  IFM  subtracted  from  the  GAIM-GM  AFWA+30  grid  specification.  The 
SED’s  location  is  clearly  indicated  by  the  dark  blue  patches  where  GAIM-GM  differed 
most  from  the  IFM.  White  regions  indicate  where  GAIM-GM  reproduced  the  IFM 
background  with  very  little  error 


the  IFM.  In  Figure  28(a)  the  white  pixels  within  the  SED  indicate  that  GAIM-GM 
fully  replicated  some  of  the  IFM  TEC  values  using  the  12  station  AFWA  grid.  The 
GAIM-GM  solution  improved  further  in  Figure  28(b)  when  30  more  stations  were 
added,  most  notably  over  the  the  central  U.S.  and  Gulf  of  Mexico. 

The  skill  scores  for  Case  2  are  given  in  Figure  29.  The  overall  trends  are  similar 
to  those  for  Case  1  but  exhibit  a  more  jagged  appearance.  Note  that  many  of  the 
sharp  drops  in  skill  score  occurred  when  the  SED  was  over  one  of  the  AFWA  grid 
stations.  At  those  times,  the  AFWA  grid  had  direct  vertical  TEC  measurements  of 
the  SED  as  it  passed  over  GPS  stations,  resulting  in  lower  errors  for  the  SED-only 
domain.  The  CORs  grids  with  additional  stations  did  not  necessarily  perform  poorly 
at  those  times,  but  rather  the  AFWA  grid  performed  exceptionally  well.  For  instance, 
whenever  the  SED  passed  over  an  AFWA  station,  the  SED  was  measured  well  and 
the  AFWA  grid  RMSEs  decreased,  which  in  turn  lowered  the  skill  scores  for  all  other 
grids  (see  Equation  8).  Also  note  in  Figure  29(c)  that  skill  scores  for  the  SED-only 
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Figure  29.  Skill  scores  for  each  of  the  CORS  grids  during  the  storm  period  from  1800 
—  2045  UT.  Compared  are  the  accuracy  improvement  over  the  AFWA  grid  for  three 
domains:  (b)  the  entire  map,  (c)  SED  only,  and  (d)  everything  outside  the  SED.  The 
trajectory  of  the  SED  is  shown  in  (a)  as  a  reference 


domain  increased  as  the  number  of  stations  increased,  whereas  in  Figure  29(d)  skill 
scores  decreased  for  everything  outside  the  SED  as  stations  increased.  This  interesting 
relationship  was  also  seen  in  Case  1. 

Table  10  summarizes  the  skill  scores  for  all  five  runs  by  averaging  the  time- 
dependent  scores  over  the  1800  -  2045  UT  storm  period.  For  the  domains  of  the 
entire  map  and  SED  only,  these  scores  were  an  average  of  10%  and  21%  higher,  re¬ 
spectively,  than  the  scores  in  Case  1.  For  the  domain  of  outside  the  SED,  they  were 
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an  average  of  15%  lower  than  Case  1.  Again  it  appeared  that  the  entire  map  domain 
saw  little  benefit  from  adding  more  than  30  stations  to  the  AFWA  grid,  whereas  the 
SED-only  domain  saw  a  steady  increase  in  accuracy.  All  areas  outside  the  SED  had  a 
steady  decline  in  accuracy  as  stations  increased.  Note  that  the  zero  station  solutions 
in  all  three  domains  were  50  -  60  skill  score  points  below  the  AFWA  grid  (compared 
to  10  -  40  points  below  the  AFWA  grid  for  Case  1).  This  demonstrated  the  modified 
version  of  GAIM-GM’s  impressive  ability  to  reproduce  the  SED  with  smaller  amounts 
of  assimilated  data. 


Table  10.  Summary  of  Case  2  average  skill  scores.  Standard  deviations  are  given  in 
parenthesis  to  the  right  of  each  score 


Domain 

No  stations 

AFWA 

+30 

AFWA 

+100 

AFWA 

+200 

AFWA 

+400 

Entire  map 
SED  only 
Outside  SED 

-60.0  (19.1) 
-64.2  (31.5) 
-56.6  (4L0) 

13.3  (6.2) 
13.1  (11.5) 
17.7  (8.3) 

15.4  (7.6) 
21.2  (18.3) 
15.4  (7.3) 

15.6  (8.0) 

26.4  (19.9) 

11.5  (7.0) 

15.4  (7.8) 
31.0  (21.7) 
7.3  (7.9) 

4.5  Case  3:  Modified  GAIM-GM  and  Ideal  Grids 

This  section  introduces  results  from  the  Ideal  grids,  comprising  21,  90,  360,  and 
720  evenly  spaced,  artificial  GPS  ground  stations  (see  Figure  18).  These  grids  were 
created  with  the  goal  of  finding  the  best  possible  GPS  receiver  distribution  over  the 
CONUS  area  without  the  location  restrictions  of  the  current  CORS  network. 

Figure  30  shows  the  TEC  results  from  the  four  GAIM-GM  model  runs  at  1945  UT, 
compared  to  the  IFM  background  with  imposed  SED.  At  this  time  the  highest  TEC 
value  in  the  IFM  synthetic  SED  was  98.9  TEC  units.  The  Ideal  21  grid  reached  a  max 
of  75.7  TEC  units  within  the  SED,  followed  by  80.1,  87.6,  and  88.1  TEC  units  for 
the  Ideal  90,  Ideal  360,  and  Ideal  720  grids,  respectively.  These  values  indicate  that 
GAIM-GM’s  accuracy  in  the  SED-only  domain  improved  as  station  count  increased, 
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as  expected.  Nevertheless  it  was  surprising  that  doubling  the  number  of  stations  from 
360  to  720  only  yielded  a  0.5  TEC  unit  increase  in  the  maximum  reproduced  TEC 
value. 

The  greatest  advantage  the  Ideal  grids  had  over  the  CORS  grid  was  the  reduction 
of  data  sparse  areas.  This  is  illustrated  in  Figure  31,  which  gives  differences  between 
the  GAIM-GM  specifications  and  the  IFM  background  for  the  AFWA  grid  (  31(a)) 
and  Ideal  90  grid  (  31(b)).  Note  particularly  the  Gulf  of  Mexico  and  Atlantic  Ocean 
areas.  The  Ideal  90  grid  had  GPS  receivers  over  these  regions  and  was  able  to  dras¬ 
tically  improve  its  specification  of  the  IFM  background,  as  opposed  to  the  CORS 
AFWA  grid. 

The  skill  score  graphs  against  the  AFWA  grid  are  given  in  Figure  32.  There  are 
a  three  key  features.  First,  since  the  Ideal  grids  covered  the  entire  GAIM-GM  model 
grid,  their  skill  scores  were  not  as  dependent  upon  the  SED’s  location  as  the  scores 
for  Cases  1  and  2  were.  The  Ideal  21  grid  was  an  exception,  since  GPS  stations 
were  so  widely  separated  that  the  SED’s  path  over  those  stations  had  a  strong  effect 
on  skill  scores.  More  on  that  topic  will  be  presented  in  Case  7.  Second,  the  skill 
scores  changed  more  drastically  between  the  two  least-dense  Ideal  grids  than  they  did 
between  the  two  least-dense  CORS  grids.  This  was  likely  due  to  the  rapid  decrease 
in  average  distance  between  stations  (as  shown  in  Table  6).  Third,  the  Ideal  360  and 
Ideal  720  grids  gave  very  similar  results  for  the  entire  map  and  SED-only  domains,  in 
spite  of  the  2x  difference  in  ground  station  density.  This  suggests  that  slant  TEC  data 
saturation  in  GAIM-GM  may  have  occurred  before  each  model  grid  square  contained 
an  observation. 

Table  11  lists  the  time-dependent  skill  scores  averaged  over  the  1800  -  2045  UT 
storm  period  for  all  four  Ideal  grids.  Owing  to  data  sparsity,  the  Ideal  21  grid  per¬ 
formed  the  worst  with  a  2.1  score  for  the  entire  map  and  -3.9  for  the  SED-only  domain. 
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Figure  30.  GAIM-GM  TEC  reproduction  of  the  synthetic  SED  at  1945  UT,  using  GPS 
slant  TEC  data  from  the  Ideal  grids.  The  synthetic  SED  in  the  IFM  background  is 
shown  in  (a) 
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Figure  31.  TEC  difference  between  GAIM-GM  and  IFM  background  at  1945  UT. 
The  IFM  was  subtracted  from  the  GAIM-GM  specification  for  the  AFWA  grid  (a) 
and  Ideal  90  grid  (b).  The  SED’s  location  is  clearly  indicated  by  the  dark  blue  patch 
where  GAIM-GM  differs  most  from  the  IFM.  White  regions  indicate  where  GAIM-GM 
reproduced  the  IFM  background  with  very  little  error 


The  Ideal  360  grid  scored  the  highest  for  the  entire  map  with  31.1  points.  Skill  scores 
for  the  SED-only  domain  steadily  increased  as  GPS  station  counts  increased,  reach¬ 
ing  a  high  of  47.0  with  the  Ideal  720  grid.  However,  the  entire  map  domain  had  a 
lower  score  with  720  stations  than  with  only  360  or  90  stations,  indicating  that  errors 
outside  the  SED  accumulated  as  station  count  increased. 


Table  11.  Summary  of  Case  3  average  skill  scores.  Standard  deviations  are  given  in 
parenthesis  to  the  right  of  each  score 


Domain 

No  stations 

Ideal  21 

Ideal  90 

Ideal  360 

Ideal  720 

Entire  map 
SED  only 
Outside  SED 

-60.0  (19.1) 
-64.2  (31.5) 
-56.6  (41-0) 

2.1  (15.5) 
-3.9  (25.4) 
8.7  (11.8) 

27.7  (7.0) 
23.9  (12.1) 

31.8  (8.8) 

31.1  (5.8) 
42.8  (6.8) 

18.1  (16.8) 

27.0  (8.9) 
47.0  (6.9) 

6.4  (23.2) 

4.6  Case  4:  High- Altitude  SED 


Cases  4  through  7  are  one-run  scenarios  (see  Table  8  for  descriptions) ,  unlike  the 
previous  cases  which  were  sets  of  runs  using  4-5  different  GPS  station  distributions. 
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Figure  32.  Skill  scores  for  each  of  the  Ideal  grids  during  the  storm  period  from  1800 
—  2045  UT.  Compared  are  the  accuracy  improvement  over  the  AFWA  grid  for  three 
domains:  (b)  the  entire  map,  (c)  SED  only,  and  (d)  everything  outside  the  SED.  The 
trajectory  of  the  SED  is  shown  in  (a)  as  a  reference 
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Each  of  these  next  cases  has  its  own  significant  variation  from  the  3  previous  cases 
and  will  shed  additional  light  on  GAIM-GM’s  abilities  to  create  accurate  TEC  maps 
during  geomagnetic  storming,  beginning  with  Case  4. 

The  motivation  for  this  case  is  illustrated  in  Figure  33,  which  shows  a  sequence 
of  GAIM-GM  vertical  density  profiles  from  Case  1  for  1945  UT.  Figure  33(a)  is  a 
specification  with  zero  data  assimilation,  Figure  33(b)  is  a  specification  with  the 
AFWA+100  grid,  and  Figure  33(c)  gives  the  difference  between  the  two.  GAIM-GM 
correctly  recreated  the  plasma  enhancement  in  the  SED  layer  from  200  -  600  km, 
but  inaccurately  increased  electron  densities  in  the  rest  of  the  vertical  column  when 
it  should  have  left  them  alone.  In  other  words,  it  put  the  SED  exactly  where  the 
F2  peak  is  climatologically  expected  to  be.  This  raised  the  question:  how  well  can 
GAIM-GM  reproduce  the  vertical  profile  of  a  SED  if  the  enhanced  densities  do  not 
fall  within  the  climatological  boundaries  of  the  F2  peak? 

To  answer  that  question,  a  special  high-altitude  SED  was  designed.  It  used  the 
same  IFM  background  but  the  SED  layer  was  set  at  500  -  1000  km  and  all  densities 
within  this  layer  were  enhanced  by  a  factor  of  4.  Densities  in  the  layer  from  90  - 
500  km  were  depleted  by  a  factor  of  2  to  simulate  the  base  of  the  SED  being  lifted 
(Figure  34(a)).  Slant  TEC  measurements  were  taken  using  the  AFWA+100  grid  and 
loaded  into  the  modified  version  of  GAIM-GM,  which  was  done  to  give  the  model  the 
best  chance  of  accurately  reproducing  the  SED. 

The  results  are  shown  in  Figure  34(b).  GAIM-GM  only  lifted  the  base  of  the  SED 
about  50  km  above  the  background,  specifying  it  near  350  km  instead  of  at  500  km.  In 
the  SED  core  GAIM-GM  reached  a  maximum  density  of  1.09xl06  electrons  /  cm~3 , 
compared  to  the  maximum  of  3.34xl06  electrons  /  cm~3  in  the  IFM. 

This  highlights  one  of  the  shortcomings  of  the  Gauss-Mar kov  version  of  GAIM, 
which  presently  does  not  include  the  horizontal  winds  nor  low-latitude  electric  fields 
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(a)  Zero  Stations 


(b)  AFWA+100 


(c)  Density  difference  (b  minus  a) 


Figure  33.  Motivation  for  Case  4,  taken  from  Case  2  vertical  density  profiles  at  40° 
N  latitude.  Shown  are  1945  UT  specifications  for:  (a)  GAIM-GM  vertical  profile 
with  zero  stations,  (b)  GAIM-GM  vertical  profile  with  AFWA+100  stations,  and  (c) 
the  difference  between  the  two.  Note  that  GAIM-GM  increased  electron  densities 
throughout  the  entire  vertical  column  and  not  just  in  the  200  —  600  km  SED  layer 
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Figure  34.  GAIM-GM  electron  density  reproduction  of  the  high  altitude  synthetic 
SED,  using  the  AFWA+100  grid  at  1945  UT.  Both  latitude  slices  were  taken  at  40°  N. 
Note  how  GAIM-GM  did  not  accurately  capture  the  vertical  profile  of  the  SED 


necessary  to  capture  the  vertical  profile  of  the  SED  ( Scherliess ,  2011).  A  physics- 
based  version  of  GAIM  is  currently  under  development  at  USU  and  is  expected  to 
significantly  improve  the  3-D  plasma  density  reconstruction  {Scherliess ,  2004). 

4.7  Case  5:  Reduced  Time  Constant 

The  Time  Constant  determines  how  long  ingested  data  are  retained  in  GAIM-GM 
as  density  perturbations.  The  default  Time  Constant  is  5  and  causes  perturbations 
from  assimilated  data  to  be  held  in  GAIM-GM’s  density  perturbation  matrix  for  5 
hours  before  being  relaxed  to  zero.  The  purpose  of  the  Time  Constant  is  to  retain 
ionospheric  features  detected  by  observations  and  propagate  them  over  data  sparse 
areas.  For  this  case  the  modified  version  of  GAIM-GM  was  run  with  the  AFWA+400 
grid  and  the  Time  Constant  was  reduced  to  1.  The  results  were  then  compared  to 
those  from  Case  2. 

Reducing  the  Time  Constant  caused  the  greatest  effects  off  the  western  coast  of 
CONUS  and  Mexico.  Figure  35  is  a  difference  plot  at  1945  UT  between  Case  2-E 
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and  Case  5  TEC  specifications,  where  Case  2-E  was  subtracted  from  Case  5.  The 
blue  shading  west  of  the  Pacific  Ocean  coastline  indicates  that  there  was  4  -  7%  less 
TEC  in  this  area  when  the  Time  Constant  was  reduced.  This  is  because  the  higher- 
value  perturbations  over  CONUS,  introduced  from  slant  TEC  data,  were  relaxed 
from  the  density  perturbation  matrix  before  they  could  propagate  very  far  from  their 
originating  GPS  ground  stations.  Smaller  TEC  differences  also  occurred  over  the 
mainland  in  the  lee  of  dense  GPS  ground  station  areas. 
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Figure  35.  TEC  difference  at  1945  UT  caused  by  lowering  the  Time  Constant  from  5  to 
1  with  the  AFWA+400  grid  in  the  modified  version  of  GAIM-GM.  The  blue,  negative 
difference  areas  occurred  where  plasma  perturbations  moved  from  areas  of  dense  GPS 
station  coverage  to  areas  of  sparse  (or  no)  GPS  station  coverage  and  relaxed  out  of  the 
model  solution  before  drifting  very  far  downstream 


Removing  these  density  perturbations  from  the  model  after  just  1  hour  instead 
of  5  gave  a  significant  boost  to  the  entire  grid’s  skill  score.  The  average  skill  score 
for  the  entire  map  over  the  1800  -  2145  UT  storm  period  increased  by  9.3  points, 
while  the  SED-only  domain  saw  an  increase  of  4.4  points.  Figure  36  shows  the  time- 


65 


(a)  SED  movement 


UT  Hour 

(b)  Skill  Score 


Figure  36.  Skill  scores  for  the  entire  map,  using  the  AFWA+400  grid  in  the  modified 
version  of  GAIM-GM.  A  trajectory  plot  of  the  synthetic  SED’s  movement  is  shown  in 
(a).  The  sharp  increase  in  scores  at  1915  UT  occurred  as  the  SED  moved  onshore  and 
over  more  ground  receivers 


dependent  skill  score  improvement  for  the  entire  grid  alongside  a  diagram  of  the  SED’s 
trajectory.  The  difference  between  scores  is  slightly  larger  at  the  end  of  the  storm 
period,  most  likely  due  to  less  high-density  SED  perturbations  advecting  downstream 
into  GPS  station-sparse  areas.  A  summary  of  average  skill  scores  for  all  three  domains 
is  presented  in  Table  12. 


Table  12.  Summary  of  Case  5  average  skill  scores.  Standard  deviations  are  given  in 
parenthesis  to  the  right  of  each  score 


Domain 

Time  Constant  =  5 

Time  Constant  —  1 

Entire  map 
SED  only 
Outside  SED 

15.4  (7.8) 

31.0  (21.7) 

7.7  (7.9) 

24.7  (8.1) 

35.4  (17.5) 

18.9  (10.2) 

4.8  Case  6:  Simulated  Bad  Data  at  One  Station 


Thus  far  the  modified  version  of  GAIM-GM  (which  trusts  data  more)  has  shown 
overall  higher  skill  scores  than  normal  version  of  GAIM-GM.  However,  all  slant  TEC 
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measurements  used  so  far  have  been  carefully  controlled  and  relatively  error-free  since 
they  are  artificial  data.  The  objective  of  this  case  was  to  deliberately  introduce 
erroneous  data  into  GAIM-GM  and  determine  how  it  affected  the  model’s  accuracy. 

The  AFWA+100  grid  was  used  for  this  case.  Freeport,  Ohio,  a  GPS  ground 
station  which  resides  directly  in  the  SED’s  path,  was  selected  to  give  bad  data  to 
GAIM-GM  for  day  172.  Every  slant  TEC  measurement  was  deliberately  reduced  by 
a  factor  of  4  to  simulate  corrupt  data,  which  in  reality  could  result  from  things  such  as 
uncalibrated  instrumentation,  broken  sensors,  or  interference  from  nearby  microwave 
radiation  sources,  among  other  things. 

Figure  37  shows  the  percentage  TEC  difference  at  1945  UT  caused  by  introducing 
the  bad  data.  This  figure  was  created  by  subtracting  the  AFWA+100  grid  output 
with  good  data  (see  Case  2-C)  from  the  Case  6  output  which  include  the  bad  station. 
The  bad  data  caused  TEC  depletions  as  high  as  9%  over  the  Appalachian  Mountains 
and  into  the  Tennessee  River  Valley,  and  a  TEC  increase  of  up  to  10%  over  the  lower 
Great  Lakes.  Note  that  the  effects  of  the  bad  data  were  not  mitigated  by  the  dozen 
other  stations  in  the  region  that  provided  good  data. 

Examining  the  entire  grid’s  skill  score  over  the  course  of  the  day  reveals  interesting 
behavior.  Note  in  Figure  38  that  bad  data  at  Freeport  caused  three  negative  skill 
score  spikes  at  1345  UT,  1745  UT,  and  2145  UT.  The  1345  UT  skill  score  decrease  was 
also  seen  to  a  lesser  extent  in  the  good  data  model  run  (dashed  line)  and  coincided 
with  a  large  TEC  gradient  from  the  Appleton  Anomaly  moving  over  the  CONUS 
East  Coast.  All  three  negative  skill  score  spikes  due  to  bad  data  (solid  line)  were 
unexpectedly  greater  than  25  points. 

The  skill  score  maps  in  Figure  39  shed  light  on  the  cause  of  one  of  the  skill  score 
decrease.  On  the  left  is  the  1730  UT  skill  score  map  and  on  the  right  is  the  1745 
UT  map.  A  large  area  of  negative  skill  scores  appeared  over  the  lower  right  corner  of 
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Figure  37.  TEC  difference  in  the  modified  version  of  GAIM-GM  at  1945  UT  caused  by 
bad  data  from  one  GPS  ground  station.  Small  circles  represent  111  stations  with  good 
data  while  the  large  black  circle  over  Ohio  is  a  station  with  deliberately  degraded  slant 
TEC  data  (reduced  by  a  factor  of  4).  Note  the  far-reaching  effects  despite  the  dozens 
of  good  data  sources  in  the  region.  The  Time  Constant  was  5 


UT  Hour 


Figure  38.  Skill  score  for  the  entire  map,  comparing  the  AFWA+100  grid  with  good 
data  to  the  AFWA+100  grid  where  one  station’s  data  were  reduced  by  a  factor  of  four 
to  simulate  erroneous  data.  The  modified  version  of  GAIM-GM  was  used.  The  storm 
period  is  defined  by  the  two  vertical  asymptotes.  Note  the  three  sudden  drops  in  skill 
score  in  the  latter  half  of  the  day  due  to  ingestion  of  bad  data 
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Figure  39.  Time-dependent  skill  score  maps  showing  the  sudden  appearance  of  a  large 
area  of  negative  skill  scores  over  the  Atlantic  Ocean.  The  modified  version  of  GAIM- 
GM  was  run  using  the  AFWA+100  list  where  data  from  one  station  over  Ohio  was 
purposely  degraded  (slant  TEC  reduced  by  a  factor  of  4).  The  area  of  abruptly  lower 
skill  scores  was  caused  by  degraded  low-angle  slant  TEC  measurements  as  two  satellites 
over  the  mid- Atlantic  Ocean  came  the  Ohio  station’s  view,  indicated  by  the  black  lines 


the  map  at  1745  UT  as  two  GPS  satellites  over  the  mid- Atlantic  Ocean  moved  into 
Freeport’s  view,  providing  two  new,  degraded  slant  TEC  measurements.  This  scenario 
could  happen  dozens  of  times  per  day  in  reality,  which  is  why  AFWA  currently  does 
not  use  this  modified  version  of  GAIM-GM  that  blindly  trusts  data.  The  normal 
version  of  GAIM-GM  is  designed  to  reject  obviously  wrong  observations  ( Scherliess 
et  al. ,  2006)  and  would  not  have  assimilated  the  bad  data  from  Freeport. 

Table  13  gives  the  average  skill  scores  during  the  1800  -  2145  UT  storm  period 
for  all  three  map  domains.  Bad  data  from  FREO  caused  the  entire  map  to  lose  1.9 
points,  the  SED  region  to  lose  2.4  points,  and  everything  outside  the  SED  to  lose  2.8 
points.  It  was  interesting  that  bad  data  from  1%  of  the  stations  could  reduce  the 
entire  map’s  skill  score  by  12%  when  modified  GAIM-GM  was  used. 


69 


Table  13.  Summary  of  Case  6  average  skill  scores.  Standard  deviations  are  given  in 
parenthesis  to  the  right  of  each  score 


Domain 

All  good  stations 

One  bad  station 

Entire  map 
SED  only 
Outside  SED 

15.4  (7.6) 

21.2  (18.3) 

15.4  (7.3) 

13.5  (7.5) 

18.8  (18-4) 

12.6  (9.2) 

4.9  Case  7:  Reversed  Ideal  21  Grid 

This  case  examined  how  the  specific  placement,  not  quantity,  of  GPS  ground 
stations  affected  GAIM-GM  TEC  maps.  Two  artificial  station  grids  were  compared: 
1)  the  Ideal  21  grid  and  2)  a  left-to-right  reversed  version  of  the  Ideal  21  grid,  hereafter 
called  the  Reversed  Ideal  21  grid. 

A  visual  comparison  of  TEC  results  at  1945  UT  for  each  grid  is  presented  in 
Figure  40,  where  Figure  40(a)  used  the  Ideal  21  grid  and  Figure  40(b)  used  the 
Reversed  Ideal  21  grid.  The  SED  was  resolved  much  better  by  the  reversed  Ideal  21 
grid  because  of  the  SED’s  specific  path  in  relation  to  the  GPS  receivers.  The  Reversed 
Ideal  21  grid  attained  a  maximum  TEC  value  in  the  SED  region  of  89.3  TEC  units 
at  1945  UT,  compared  to  only  55.8  TEC  units  attained  by  the  Ideal  21  grid  (note 
that  the  synthetic  SED  maximum  was  111.4  TEC  units). 

Figure  41(a)  shows  the  exact  path  of  the  SED  through  both  grids,  where  the  small 
black  circles  represent  Ideal  21  grid  stations  and  the  large  red  circles  are  Reversed 
Ideal  21  grid  stations.  To  its  right  in  Figure  41(b)  are  the  time-dependent  skill  scores 
for  the  entire  map.  Scores  significantly  improved  with  the  reversed  Ideal  21  grid, 
especially  prior  to  1930  UT  when  the  SED  was  still  over  the  ocean. 

The  average  skill  scores  for  all  domains  are  given  in  Table  14.  All  areas  saw 
improvements,  with  the  largest  gains  (16.5  point  improvement)  occurring  for  the  SED- 
only  domain.  This  case  demonstrated  that  21  well-placed  land  and  sea  GPS  receivers 
had  a  13.5  point  skill  score  improvement  over  the  12  AFWA  ground  receivers. 


70 


TEC  Units  TEC  Units 


(a)  Ideal  21  grid  (b)  Reversed  Ideal  21  grid 

Figure  40.  GAIM-GM  TEC  maps  at  1945  UT  using  (a)  the  Ideal  21  grid  and  (b)  a 
left-right  reversed  Ideal  21  grid.  The  modified  version  of  GAIM-GM  was  used.  By 
simply  adjusting  the  arrangement  of  GPS  ground  stations  (and  without  changing  the 
average  distance  between  stations),  the  SED  was  reproduced  more  accurately 
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(b)  Skill  Scores 


Figure  41.  Skill  scores  for  the  entire  map,  comparing  the  Ideal  21  grid  to  the  Reversed 
Ideal  21  grid.  As  a  reference,  both  grids  are  superimposed  on  the  SED  movement  plot 
in  (a).  The  small  black  circles  denote  the  Ideal  21  grid  and  the  large  red  circles  denote 
the  Reversed  Ideal  21  grid.  The  Reversed  grid’s  station  locations  were  more  favorable 
for  this  SED’s  specific  trajectory,  resulting  in  largely  improved  skill  scores 
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Table  14.  Summary  of  Case  7  average  skill  scores.  Standard  deviations  are  given  in 
parenthesis  to  the  right  of  each  score 


Domain 

Ideal  21  grid 

Reversed  Ideal  21  grid 

Entire  map 
SED  only 
Outside  SED 

2.1  (15.5) 

-3.9  (25.4) 

8.7  (11.8) 

13.5  (9.2) 

12.6  (15.0) 

11.5  (12.0) 

4.10  Case  8:  Normal  GAIM-GM  with  Depleted  IFM 

Geomagnetic  storms  often  cause  TEC  increases  such  as  SEDs,  but  they  can  also 
cause  TEC  depletions  in  other  areas  outside  SEDs  due  to  complex  electric  held  and 
upper  atmosphere  dynamics  ( Foster ,  1993).  The  objective  for  this  case  was  to  dupli¬ 
cate  the  storm  shown  in  Figure  12  by  depleting  IFM  background  densities  in  areas 
outside  the  SED.  For  example,  note  the  TEC  depletions  east  and  west  of  the  historic 
SEDs  in  Figures  12  and  13. 

Using  the  historic  SED  specifications  as  a  pattern,  all  densities  above  30°  N  lat¬ 
itude  in  the  IFM  background  were  depleted  by  a  factor  of  4,  except  for  those  in  the 
synthetic  SED.  Slant  TEC  measurements  were  then  taken  in  the  normal  fashion  using 
the  CORS  grids  and  ingested  by  normal  GAIM-GM.  An  example  of  the  depleted  IFM 
background  TEC  with  the  imposed  SED  at  1945  UT  is  shown  in  Figure  42(a),  while 
Figures  42(b)  through  42(f)  give  the  GAIM-GM  TEC  specifications  using  the  five 
CORS  grids. 

Note  in  Figure  42  that  the  TEC  across  the  entire  map  changed  very  little  once 
the  amount  of  ground  stations  exceeded  the  AFWA+30  point.  Additionally,  note 
the  green-colored  ridge  of  higher  TEC  values  along  the  eastern  edge  of  the  map 
where  GAIM-GM  lacked  the  data  to  reproduce  the  lower  background  densities.  This 
ridge  can  be  seen  in  Figure  43,  which  shows  1945  UT  electron  density  vertical  cross- 
sections  at  40°  N  latitude  for  the  synthetic  SED  in  IFM,  compared  to  the  GAIM-GM 
specifications  using  the  AFWA  and  AFWA+30  grids. 
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Figure  42.  GAIM-GM  TEC  reproduction  of  the  synthetic  SED  at  1945  UT,  using 
GPS  slant  TEC  data  from  the  CORS  grids.  The  synthetic  SED  in  the  depleted  IFM 
background  is  shown  in  (a)  as  a  reference  and  the  12  AFWA  ground  stations  are  plotted 
as  black  circles  in  (b) 


73 


1200 


(a)  IFM 


0  0.5  1  1.5  2  0  0.5  1  1.5  2 

Electron  Density  xio6  Electron  Density  xio6 


(b)  GAIM-GM:  AFWA  grid 


(c)  GAIM-GM:  AFWA+30  grid 


Figure  43.  IFM  and  GAIM-GM  electron  density  vertical  profiles  at  1945  UT.  Latitude 
slices  at  40°  N  are  shown  for  the  IFM  background  (a),  AFWA  grid  output  (b)  and 
AFWA+30  grid  output  (c) 


In  the  IFM,  the  maximum  density  of  the  SED  core  was  2.29xl06  electrons /cm~3 , 
whereas  for  the  AFWA  and  AFWA+30  grids  it  was  F56xl06  electrons  /  cm~3  and 
1.90  xlO6  electrons / cm~3 ,  respectively.  Electron  density  profiles  for  the  AFWA+100, 
AFWA+200,  and  AFWA+400  grids  were  similar  to  Figure  43(c)  and  are  not  shown. 

Skill  scores  for  the  SED-only  domain,  along  with  a  trajectory  plot  of  the  SED,  are 
shown  in  Figure  44.  Unlike  the  relatively  smooth  curves  of  Figure  26(c)  from  Case  1, 
which  did  not  have  a  depleted  IFM,  these  scores  fluctuated  sharply  with  each  time 
step.  This  pattern  was  explained  in  Case  2  as  resulting  from  very  low  RMSEs  for  the 
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AFWA  grid  whenever  the  SED  was  directly  over  an  AFWA  station.  As  expected,  the 
skill  scores  improved  as  station  counts  increased. 


Figure  44.  SED  skill  scores  over  the  entire  map  for  each  of  the  CORS  grids  during  the 
storm  period  from  1800  —  2045  UT.  The  trajectory  of  the  SED  is  shown  in  (a)  as  a 
reference.  Scores  above  0  indicate  the  grid  was  more  accurate  at  reproducing  the  IFM 
background  than  the  AFWA  grid,  while  scores  below  0  indicate  worse  accuracy 


Table  15  lists  the  average  skill  scores  during  the  1800  -  2045  UT  storm  period 
for  the  three  domains.  The  zero  station  scores  were  extremely  low  due  to  the  large 
difference  between  the  depleted  IFM  background  and  GAIM-GM’s  internal  IFM  back¬ 
ground.  This  re-emphasizes  the  advantage  that  GAIM-GM  has  over  the  IFM  during 
geomagnetic  storming  due  to  data  assimilation.  Skill  scores  increased  in  every  domain 
as  more  stations  were  added,  with  the  greatest  improvements  seen  in  the  SED-only 
domain.  With  respect  to  the  entire  map,  going  from  30  extra  stations  to  400  extra 
stations  only  caused  a  2.7  point  skill  score  increase. 
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Table  15.  Summary  of  Case  8  average  skill  scores.  Standard  deviations  are  given  in 
parenthesis  to  the  right  of  each  score 


Domain 

No  stations 

AFWA 

+30 

AFWA 

+100 

AFWA 

+200 

AFWA 

+400 

Entire  map 
SED  only 
Outside  SED 

-77.4  (14.7) 
-9.8  (18.3) 
-81.8  (46.5) 

7.1  (2.1) 

11.6  (12.7) 

7.2  (9.6) 

7.8  (2.1) 

16.8  (17.1) 
7.7  (8.6) 

8.1  (2.3) 

20.8  (16.8) 
7.8  (8.0) 

9.9  (2.3) 

24.1  (17.1) 
9.4  (8.8) 

4.11  Case  9:  Reduced  Time  Constant  with  Depleted  IFM 

This  case  tested  the  effects  of  reducing  the  Time  Constant  in  normal  GAIM-GM 
with  a  depleted  IFM  background  for  all  densities  above  30°  N  latitude.  As  mentioned 
in  Case  5,  the  purpose  of  the  Time  Constant  is  to  retain  ionospheric  features  detected 
by  observations  and  propagate  them  over  data  sparse  areas.  Recall  that  Case  5  showed 
how  reducing  the  Time  Constant  to  1  raised  skill  scores  (i.e.  reduced  the  errors)  for 
all  areas  of  the  map.  Skill  scores  were  not  available  for  this  case  because  comparing 
the  AFWA  grid  to  itself  results  in  a  skill  score  of  zero  (see  Equation  8),  so  MAEs 
were  used  instead. 

The  TEC  difference  caused  by  lowering  the  Time  Constant  is  shown  in  Figure  45. 
Note  that  this  difference  was  much  greater  than  that  seen  in  Figure  35  from  Case  5. 
Because  the  slant  TEC  measurements  through  the  depleted  IFM  background  were 
so  different  from  GAIM-GM’s  internal  IFM  background,  the  resulting  perturbations 
were  also  much  larger.  Removing  those  perturbations  after  just  1  hour  instead  of  5 
naturally  caused  the  large  TEC  differences.  Another  item  to  note  is  that  the  areas 
of  positive  and  negative  differences  nearly  balanced  each  other  out,  leading  to  the 
small  average  MAE  differences  listed  in  Table  16.  Overall  there  was  very  little  change 
in  the  storm  period’s  average  MAE  other  than  a  2.8  TEC  unit  improvement  for  the 
SED-only  domain. 


76 


<D 

-o 

D 


CO 


230  240  250  260  270  280  290  300  310 


Longitude 


-100  -80  -60  -40  -20  0  20  40  60  80  100 

Percent 


Figure  45.  TEC  difference  at  1945  UT  caused  by  lowering  the  Time  Constant  from  5 
to  1  for  the  AFWA  grid.  All  background  IFM  densities  above  30°  were  depleted  by 
a  factor  of  4  prior  to  slant  TEC  piercing,  therefore  these  large  percentage  differences 
equate  to  relatively  small  TEC  unit  differences.  Note  that  most  of  the  differences 
occured  downstream  of  the  GPS  ground  stations 


Table  16.  Summary  of  Case  9  average  MAEs  (TEC  units).  Standard  deviations  are 
given  in  parenthesis  to  the  right  of  each  score 


Domain 

Time  Constant  =  5 

Time  Constant  =  1 

Entire  map 
SED  only 
Outside  SED 

12.3  (1.9) 

15.4  (10.4) 

11.9  (3.6) 

12.3  (1.7) 

12.6  (9.6) 

11.9  (3.6) 

4.12  Case  10:  Elevation  Mask  Raised  to  45° 

The  final  case  in  this  project  examined  how  changing  the  GPS  station  Elevation 
Mask  affected  the  accuracy  of  GAIM-GM  TEC  specifications.  For  all  nine  previous 
cases  the  Elevation  Mask  was  set  to  15°,  meaning  slant  TEC  measurements  would 
only  be  taken  between  a  GPS  satellite  and  ground  station  if  the  satellite  was  at  least 
15°  above  the  horizon,  from  the  station’s  reference  frame.  In  this  case  the  Elevation 
Mask  was  raised  to  45°,  significantly  narrowing  the  field  of  view  for  every  station. 
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Slant  TEC  measurements  were  then  taken  using  the  AFWA  grid.  As  with  the  two 
prior  cases,  the  normal  version  of  GAIM-GM  was  used  and  the  IFM  background  was 
depleted  above  30°  N  latitude  by  a  factor  of  four.  The  Time  Constant  was  reset  to 
its  default  value  of  5. 

Figure  46  shows  the  outcome  of  this  scenario  by  comparing  its  results  to  electron 
density  and  TEC  maps  at  1945  UT  from  Case  8-A.  The  central  column  gives  the  15° 
Elevation  Mask  GAIM-GM  specifications  from  8-A  while  the  right  column  gives  the 
45°  results.  Perhaps  the  most  striking  feature  is  seen  in  Figure  46(c),  where  there 
is  a  circular  region  of  low  TEC  (about  30  TEC  units)  directly  over  the  Bermuda 
GPS  receiver.  With  higher  look  angles,  there  were  fewer  long-path  slant  TEC  mea¬ 
surements  taken  near  the  Earth’s  surface  to  help  expand  Bermuda’s  influence  on  the 
surrounding  region.  The  same  phenomenon  was  seen  over  southern  CONUS  where 
higher  TEC  ridges  appeared  to  wrap  around  the  periphery  of  GPS  ground  stations, 
creating  a  tunnel  vision  effect  due  to  the  reduced  held  of  view. 

Average  MAEs  for  this  case  are  given  in  Table  17,  which  indicates  that  the  MAE 
for  the  entire  map  only  increased  by  2.3  TEC  units,  or  approximately  19%.  This 
error  would  be  expected  to  decrease  if  a  higher  density  GPS  station  grid  were  used 
and  the  fields  of  view  were  able  to  overlap. 

Table  17.  Summary  of  Case  10  average  MAE  (TEC  units).  Standard  deviations  are 
given  in  parenthesis  to  the  right  of  each  score 


Domain 

Look  angle  =15° 

Look  angle  =  45° 

Entire  map 
SED  only 
Outside  SED 

12.3  (1.9) 

26.6  (10.4) 

11.9  (3.6) 

14.6  (2.8) 

29.3  (11.9) 

14.1  (4.4) 

Table  18  compares  the  number  of  data  ingested  by  GAIM-GM  for  45°  and  15° 
Elevation  Masks.  When  the  Elevation  Mask  was  raised  to  45°  the  average  amount  of 
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the  IFM  backgrounds  that  GAIM-GM  was  recreating,  which  had  been  depleted  by  a  factor  of  4  above  30°  N 
latitude.  The  central  column  shows  GAIM-GM  output  using  the  default  15°  Elevation  Mask.  The  right  column 
shows  GAIM-GM  output  using  a  45°  Elevation  Mask 


ingested  data  per  time  step  decreased  by  55%.  Surprisingly  this  had  relatively  little 
effect  on  the  average  MAE  for  the  storm  period  of  1800  -  2045  UT. 


Table  18.  Number  of  slant  TEC  measurements  assimilated  by  GAIM-GM  using  the 
AFWA  grid  and  varying  Elevation  Masks,  in  the  normal  version  of  GAIM-GM 


Time  (UT) 

Elevation  Mask  =  45° 

Elevation  Mask  =15° 

1800 

45 

77 

1815 

44 

76 

1830 

40 

71 

1845 

43 

76 

1900 

36 

79 

1915 

35 

77 

1930 

42 

73 

1945 

43 

77 

2000 

45 

75 

2015 

45 

82 

2030 

44 

80 

2045 

43 

76 

80 


V.  Conclusions  and  Recommendations 


5.1  Chapter  Overview 

This  chapter  consists  of  two  sections.  The  first  section  summarizes  the  results 
from  using  GAIM-GM  to  reproduce  a  synthetic  SED  with  varying  numbers  of  GPS 
ground  stations  and  for  several  different  model  configurations.  The  second  section 
presents  recommendations  for  future  research  with  the  GAIM-GM  model. 

5.2  Conclusions 

The  accuracy  of  TEC  maps  produced  by  GAIM-GM  varied  considerably,  based  on 
both  amount  and  distribution  of  available  GPS  ground  stations.  Using  skill  scores  as 
a  basis  for  accuracy,  Figure  47  summarizes  GAIM-GM’s  performance  in  reproducing 
the  IFM  background  with  the  imposed  synthetic  SED.  Scores  are  given  for  three 
domains:  1)  entire  map;  2)  SED-only;  and  3)  area  outside  the  SED.  Note  that  these 
scores  are  the  averaged,  time-dependent  skills  scores  for  each  grid  during  the  1800  - 
2045  UT  storm  period. 

Model  runs  using  slant  TEC  measurements  from  a  network  of  CORS  (Figure  47(a)) 
showed  TEC  skill  score  improvements  of  13.1  points  for  the  entire  map  when  30  ground 
stations  were  added  to  the  existing  AFWA  network.  Increasing  the  number  of  stations 
beyond  30  provided  less  than  1  point  of  improvement  for  the  entire  map.  However, 
skill  scores  for  the  SED-only  domain  continued  to  improve  as  stations  were  added 
and  reached  a  score  of  24.1  with  400  extra  stations.  Skill  scores  for  all  areas  outside 
the  SED  simultaneous  declined  at  about  the  same  rate.  Similar  skill  score  trends  for 
the  entire  map  and  SED-only  domain  (Figure  47(d))  were  observed  when  the  IFM 
background  was  depleted  by  a  factor  of  4  above  30°  N  latitude,  although  scores  for 
areas  outside  the  SED  did  improve  2.2  points  with  400  extra  stations. 
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Figure  47.  GAIM-GM  skill  score  summaries  for  the  4  most  important  cases,  which  are: 
(a)  normal  GAIM-GM  with  CORS  grids;  (b)  modified  GAIM-GM  with  CORS  grids; 
(c)  modified  GAIM-GM  with  Ideal  grids;  and  (d)  normal  GAIM-GM  with  CORS  grids 
and  depleted  IFM  background; 


CORS  grids  in  the  modified  version  of  GAIM-GM  performed  on  average  10% 
better  for  the  entire  grid,  21%  better  for  the  SED-only  domain,  and  15%  worse  for 
everything  outside  the  SED  (see  Figure  47(b)),  compared  to  the  normal  version  of 
GAIM-GM.  Although  the  improvements  for  the  entire  grid  and  SED-only  domain 
were  impressive,  it  must  be  emphasized  that  the  data  in  this  research  were  quality- 
controlled  and  only  had  5,7%  uncertainty.  When  slant  TEC  data  at  one  GPS  ground 
station  in  the  SED’s  path  were  deliberately  degraded  (75%  error  introduced),  skill 
scores  for  the  entire  map  dropped  an  average  of  12%.  Because  of  the  negative  effects 
even  a  small  amount  of  bad  data  can  cause,  the  modified  version  of  GAIM-GM  should 
not  be  used  operationally  as  a  stand-alone  model.  However,  it  would  be  advantageous 
to  use  it  alongside  the  normal  version  of  GAIM-GM  in  ensemble  fashion,  giving  AFWA 
access  to  potentially  superior  specification  as  long  as  data  fidelity  is  monitored. 

Tests  with  the  Ideal  grids  in  modified  GAIM-GM  (Figure  47(c))  demonstrated 
that  strategically  placing  ground  GPS  receivers  yielded  significant  improvements  in 
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accuracy.  An  evenly  distributed  grid  of  90  artificial  stations  yielded  skill  scores  of  27.1 
for  the  entire  map,  23.9  for  the  SED-only  domain,  and  31.8  for  everything  outside 
the  SED.  Those  scores  are  80%,  13%,  and  106%  higher  than  their  counterpart  scores 
using  the  AFWA+100  grid  in  modified  GAIM-GM.  It  was  also  found  that  low-density 
ground  station  grids  were  highly  dependent  upon  the  trajectory  of  the  SED.  The  Ideal 
21-station  grid  initially  scored  -3.9  for  the  SED-only  domain,  but  reversing  the  grid 
from  left  to  right  (without  altering  station  spacing)  raised  the  score  to  12.6  by  placing 
more  ground  stations  in  favorable  locations  to  observe  the  SED. 

The  significant  accuracy  advantages  of  the  Ideal  grids  could  be  realized  if  ocean- 
borne  GPS  receivers  were  added  to  AFWA’s  network.  One  possible  avenue  is  to 
deploy  GPS  receivers  on  naval  vessels,  either  on  U.S. -owned  military  ships  or  on 
commercial  vessels  via  a  contractor.  Oil  rigs  are  also  potential  sites,  although  they 
would  not  provide  as  broad  a  reach  from  the  mainland.  Buoys  might  not  be  a  viable 
option  without  a  way  to  correct  for  the  constant  altitude  huctuations  that  ocean  waves 
would  cause  to  the  receiver’s  position  (recall  that  the  GPS  receiver’s  coordinates  must 
be  precisely  known). 

In  addition  to  altering  the  number  and  location  of  ground  GPS  receivers,  adjusting 
the  Time  Constant  affected  the  accuracy  of  GAIM-GM  TEC  maps.  Reducing  the 
Time  Constant  from  5  to  1  caused  GAIM-GM  to  relax  density  perturbations  sooner 
and  had  the  greatest  effects  downstream  of  data-rich  areas.  Skill  scores  in  the  modified 
version  of  GAIM-GM  increased  by  an  average  of  9.3  for  the  entire  grid,  4.4  for  the 
SED-only  domain,  and  11.2  for  everything  outside  the  SED.  In  normal  GAIM-GM 
there  were  large  localized  TEC  differences  for  the  entire  grid,  but  negligible  overall 
reduction  of  errors.  Average  MAEs  for  the  SED-only  domain  decreased  by  3.6  TEC 
units  (a  10%  improvement  in  MAE)  because  GAIM-GM  propagated  fewer  low-density 
upstream  perturbations  into  the  SED. 
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Lowering  the  Time  Constant  in  the  operational  GAIM-GM  will  yield  accuracy 
improvements  in  some  cases  but  degradations  in  others.  With  the  Time  Constant 
set  to  1,  short-lived  features  or  those  that  have  poleward  velocity  components  would 
not  be  overly  skewed  to  the  west,  and  therefore  GAIM-GM  would  have  improved 
precision.  Conversely,  long-lived  ionospheric  features  that  propagate  zonally  west¬ 
ward,  especially  over  data-sparse  regions,  would  not  be  captured  as  well  with  a  Time 
Constant  of  1.  Given  that  AFWA  currently  uses  a  low-density  GPS  station  grid  over 
CONUS,  there  would  be  little  operational  advantage  in  lowering  the  Time  Constant. 

Changing  the  Elevation  Mask  also  affected  GAIM-GM’s  performance.  Raising 
the  angle  from  the  default  15°  to  45°  with  the  AFWA  grid  effectively  narrowed  each 
GPS  ground  station’s  field  of  view  and  limited  the  amount  of  available  slant  TEC 
observations.  The  average  TEC  MAE  increased  by  16%  for  all  areas  as  a  result  of 
this  tunnel  vision  effect,  indicated  poorer  model  performance. 

One  limitation  of  GAIM-GM  uncovered  during  this  research  was  that,  although  it 
typically  reproduced  TEC  with  good  accuracy,  it  could  not  correctly  specify  the  elec¬ 
tron  density  vertical  profile  of  a  high- altitude  synthetic  SED.  It  appeared  that  GAIM- 
GM  tended  to  reconstruct  the  height  of  the  ionosphere  F2  peak  based  on  climatology. 
Because  the  test  only  used  slant  TEC  measurements,  it  is  unclear  whether  those 
data  alone  provided  inadequate  information  for  resolving  the  SED  profile,  or  if  it  was 
entirely  due  to  the  limitations  of  the  Gauss-Markov  model.  Assimilating  ionosonde 
and  in  situ  DMSP  data  into  GAIM-GM  could  potentially  supply  the  altitude-specific 
electron  density  information  needed  for  GAIM-GM  to  accurately  specify  the  SED.  It 
is  recommended  that  Case  4  be  repeated  with  these  additional  data  types. 

Based  on  all  results,  there  are  three  recommendations  for  the  number  of  ground 
GPS  receivers  AFWA  should  employ  over  CONUS  to  help  ensure  reliable  TEC  correc¬ 
tion  maps  during  geomagnetic  storming.  These  numbers  consider  the  balance  between 
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accuracy  over  the  entire  grid  and  accuracy  for  the  SED-only  domain.  For  a  network 
using  existing  CORS  sites,  the  optimal  number  is  ~110  GPS  receivers.  If  availability 
of  data  or  economic  considerations  preclude  using  that  amount,  the  number  could  be 
lowered  to  ^40  and  still  provide  significant  improvements  over  the  current  11-station 
AFWA  network.  If  the  Ideal  grid  were  attainable  and  GPS  receivers  could  be  evenly 
situated  over  land  and  sea,  the  recommended  number  would  be  ~90. 

Finally,  note  that  these  results  were  obtained  from  very  specific  sets  of  data, 
under  carefully  controlled  scenarios,  and  that  the  synthetic  observations  were  taken 
at  GAIM-GM’s  highest  allowable  cadence  of  15  minutes.  Because  real  slant  TEC 
data  are  subject  to  many  uncertainties  (instrument  calibration,  sensor  errors,  etc.) 
and  are  not  alway  available  every  15  minutes,  results  using  real  data  in  GAIM-GM 
may  vary  from  those  presented  here. 

5.3  Recommendations  for  Future  Work 

Questions  raised  during  this  research  could  provide  several  avenues  for  future 
study  on  SEDs  and  their  effect  on  geo-location  accuracy.  These  include  repeating 
some  of  the  most  interesting  cases  but  under  different  conditions,  such  as:  1)  running 
the  Ideal  grids  in  normal  GAIM-GM  instead  of  the  modified  version  of  GAIM-GM;  2) 
conducting  tests  with  the  saturated  Ideal  720  grid  to  determine  why  it  gave  skill  scores 
lower  than  the  Ideal  360  grid;  3)  testing  Time  Constant  changes  is  sparser  grids,  such 
as  the  AFWA+30  and  AFWA+100  grids;  4)  running  more  scenarios  with  simulated 
bad  data  to  determine  why  reduced  slant  TEC  values  over  Freeport  caused  counter¬ 
intuitive  TEC  increases  over  the  Great  Lakes;  and  5)  recreating  a  high-altitude  SED 
using  DMSP  and  ionosonde  data,  as  mentioned  in  the  previous  section. 

It  is  highly  recommended  that  these  suggested  studies,  or  any  similar  studies,  take 
advantage  of  historical  SED  specifications  in  GAIM-GM  and  use  an  actual  GAIM- 
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GM  density  background  instead  of  an  IFM  background  for  slant  TEC  piercing.  Recall 
that  the  synthetic  SED  in  this  project  was  basic  construct  representing  a  single  prop¬ 
agating  patch  of  enhanced  plasma,  precisely  imposed  onto  an  IFM  background.  This 
simple  yet  effective  storm  feature  could  be  made  even  more  realistic  and  dynamic  by 
adopting  the  actual  GAIM-GM  historical  recreation  of  a  SED  as  the  density  back¬ 
ground,  and  then  taking  slant  TEC  measurements  through  it  using  a  modified  version 
of  SlantTEC.m.  With  the  necessary  MATLAB®utilites  already  pioneered  by  AFIT, 
research  could  commence  almost  immediately. 

Additionally,  the  same  cases  presented  here  could  be  repeated  using  smaller  GPS 
station  increments,  such  as  10,  20,  30,  40,  and  50.  Given  the  cost  of  GPS  data  vs. 
economic  realities,  it  is  extremely  likely  that  additions  to  the  AFWA  station  network 
will  come  in  small  increments,  perhaps  even  one  station  at  a  time.  A  methodical 
study  of  SED  climatology  to  identify  regions  where  GPS  receivers  are  needed  most, 
followed  by  testing  in  GAIM-GM,  is  recommended.  This  would  prioritize  which  GPS 
ground  stations  to  acquire  next  and  would  be  tremendously  valuable  to  AFWA  and 
space  weather  forecasting  in  the  coming  years. 


Appendix  A.  The  Global  Positioning  System  (GPS) 


The  GPS  is  a  U.S.  government- built  system  consisting  of  three  segments:  1)  space 
segment;  2)  user  segment;  and  3)  control  segment.  The  space  segment  is  operated  and 
maintained  by  the  U.S.  Air  Force  and  is  composed  of  a  constellation  of  32  satellites 
orbiting  at  approximately  20,200  km  with  an  orbital  period  of  12  hours.  Each  satellite 
transmits  two  signals,  LI  and  L2,  at  microwave  frequencies.  The  LI  frequency  is 
1575.42  MHz  and  transmits  the  Standard  Positioning  Service  signal,  which  can  be 
used  worldwide  by  anyone  without  restrictions  or  charges.  The  L2  signal  is  broadcast 
encrypted  at  1227.60  MHz  and  is  used  by  the  U.S.  military  and  its  allies’  militaries, 
as  well  as  other  select  agencies  approved  by  the  U.S.  Government  ( Dana ,  2010). 

Figure  48  shows  the  orbital  tracks  of  all  32  GPS  satellites  over  a  24-hour  period. 
Note  how  the  orbital  inclination  prevents  overflight  of  the  poles  but  still  offers  ample 
coverage  for  the  mid-latitudes  where  SEDs  frequently  occur. 

The  user  segment  consists  of  the  GPS  receiver.  At  any  given  time  there  are 
five  to  eight  satellites  with  direct  lines-of- sight  to  a  user  at  any  point  on  earth,  al¬ 
though  linkage  is  only  needed  with  four  satellites  in  order  to  determine  the  user’s 
three-dimensional  geo-location.  Three  satellites  are  used  to  determine  the  latitude, 
longitude,  and  elevation  of  the  receiver  while  the  fourth  satellite  corrects  the  time 
offset  in  the  receiver’s  clock.  Time  correction  is  necessary  so  that  the  receiver  can 
precisely  measure  the  time  delays  of  the  other  three  satellites’  signals.  If  the  receiver 
equipment  contained  its  own  atomic  clock,  perfectly  synchronized  with  the  GPS  con¬ 
stellation,  it  would  only  need  three  satellites  for  a  location  fix. 

The  third  segment  is  the  control  segment,  which  composed  of  the  Master  Control 
Facility  located  at  Schriever  Air  Force  Base,  Colorado  and  a  worldwide  monitoring 
network  of  GPS  ground  tracking  stations.  The  control  segment  performs  a  number 
of  vital  functions:  1)  monitors  and  maintains  the  satellites  in  their  proper  orbits 
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Figure  48.  Orbital  tracks  for  all  32  GPS  satellites  over  the  course  of  one  day,  at 
15-minute  time  steps.  The  orbital  inclination  of  55°  prevents  overflight  of  the  polar 
regions 

through  occasional  command  maneuvers;  2)  adjusts  the  satellites’  clocks;  3)  tracks 
the  satellites;  4)  uploads  updated  navigational  data;  and  5)  maintains  the  health  and 
status  of  the  satellite  constellation.  There  are  currently  more  than  360  active  sta¬ 
tions  in  the  International  Global  Navigation  Satellite  System  Service  (IGS)  tracking 
network  ( International  Global  Navigation  Satellite  System  Service ,  2011).  Closer  to 
CONUS,  the  National  Geodetic  Survey  (NGS)  manages  a  network  of  over  1900  Con¬ 
tinuously  Operating  Reference  Stations  (CORS)  throughout  the  U.S.,  its  territories, 
and  a  few  foreign  countries  ( National  Geodetic  Survey ,  2011).  Figure  49  shows  the 
CORS  network  over  the  CONUS  area. 

In  most  cases  the  ionosphere  is  a  hindrance  to  radio  transmission  between  GPS 
satellites  and  ground  stations  ( Ondoh  and  Marubashi ,  2001).  At  the  same  time, 
although  GPS  signals  were  not  intended  specifically  for  ionospheric  monitoring,  they 
can  be  used  for  ionospheric  research.  GPS  ground  stations  contribute  to  studies  of  the 
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Figure  49.  Continuously  Operating  Reference  Stations  (CORS)  network  over  CONUS. 
There  are  1920  stations  shown  here,  with  the  majority  located  over  the  eastern  U.S. 

ionosphere  by  calculating  slant  TEC  values,  which  are  integrated  electron  densities 
along  the  path  between  GPS  satellites  and  the  ground  stations.  The  development 
of  the  calculation  presented  here  closely  mirrors  that  given  by  Ondoh  and  Marubashi 
(2001)  and  involves  measuring  the  time  delay  difference  5t  between  the  LI  and  L2 
signals  received  at  the  station: 


St  =  A  t2  —  At  i 


(10) 


The  variables  Ati  and  A t2  correspond  to  LI  and  L2  and  represent  the  amount  of 
group  delay  imposed  on  a  radio  signal  of  frequency  /i  or  f2.  Generically  this  is: 


At  = 


8TT2mee0f2 


(11) 


where  the  constants  are  listed  in  Table  19.  The  integral  represents  the  integrated 


electron  density  over  the  propagation  path  length  (from  satellite  to  ground) : 


nds  =  Nt 


(12) 


where  NT  is  the  TEC  along  the  path  length.  To  solve  for  the  TEC,  Equations  (11) 
and  (12)  are  substituted  into  Equation  (10),  resulting  in: 


5t  = 


87r2mee0 


(13) 


Finally,  the  values  of  the  constants  in  Table  19  are  inserted  into  Equation  (13): 


Nt  =  2.85  x  1025  5t  [electrons / m2 


(14) 


From  Equation  (14)  the  slant  TEC  can  be  easily  calculated  for  any  given  amount  of 
group  delay. 


Table  19.  Slant  TEC  equation  constants 


Symbol 

Name 

Value 

e 

Elementary  charge 

1.6022  x  10“19  C 

me 

Static  electron  mass 

9.1094  x  10-31  kg 

e0 

Permittivity  of  free  space 

8.8542  x  10-12  F/m 

fi 

Lx  signal  frequency 

1575.42  MHz 

h 

L2  signal  (P  code)  frequency 

1227.60  MHz 
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Appendix  B.  AFWA  12th  Station 


The  original  list  of  GPS  ground  stations  used  by  AFWA  for  running  GAIM-GM 
contained  59  global  stations,  11  of  which  were  located  in  the  CONUS  assimilation 
region  (see  Table  3).  This  list  was  provided  courtesy  of  AFWA  and  identified  all 
sites  which  had  provided  slant  TEC  measurements  to  GAIM-GM  during  the  month 
of  February  2011.  The  actual  number  of  stations  used  changes  slightly  from  month- 
to-month,  based  on  the  availability  of  data  ( Reich ,  2011).  For  example,  during  May 
2011  GAIM-GM  used  77  stations  from  around  the  globe,  which  were  18  more  than 
the  number  used  in  February. 

Because  the  AFWA  list  has  so  few  GPS  ground  stations  over  CONUS,  there  are 
data-sparse  areas  which  prevent  accurate  specifications  of  the  ionosphere  via  slant 
TEC  measurements.  Large  scale  features  might  be  big  enough  to  blanket  a  few  GPS 
ground  stations,  but  small  scale  features  may  be  partially  or  fully  unresolved  if  they 
fall  inconveniently  between  widely  separated  stations.  Conversely,  a  small  feature 
might  pass  directly  over  an  isolated  GPS  ground  station  and  cause  GAIM-GM  to 
correlate  the  feature’s  density  into  the  surrounding  grid  squares’  density  solution, 
effectively  smearing  the  feature  over  an  extended  area  where  it  should  not  exist. 

Consider  the  two  TEC  maps  in  Figure  50.  On  the  left  is  a  GAIM-GM  historical 
specification  of  the  20  November  2003  SED.  On  the  right  is  an  IFM  specification  with 
an  imposed  synthetic  SED.  The  large  surface  area  of  enhanced  TEC  in  the  real  SED 
covers  nearly  the  bottom  third  of  the  map.  Large  Latitude  Correlations  and  Time 
Constants  would  favor  the  real  SED  in  a  data  sparse  environment  (like  the  AFWA 
grid)  by  spreading  enhanced  densities  over  areas  lacking  ground  stations,  creating 
a  more  accurate  specification.  With  the  synthetic  SED,  the  aim  is  to  minimize  the 
amount  of  smearing  caused  by  the  Latitude  Correlation  and  Time  Constant,  so  the 
best  path  is  one  through  a  corridor  of  stations  that  will  resolve  all  edges  of  the  SED. 
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Figure  50.  Historical  depiction  of  an  SED  in  GAIM-GM  (left),  versus  a  synthetic  SED 
(right).  The  synthetic  SED’s  small  profile  make  it  a  challenging  target  for  slant  TEC 
measurements  over  sparsely  populated  GPS  ground  station  networks 


In  the  case  of  the  original  11-station  AFWA  grid,  the  synthetic  SED  happens  to 
pass  directly  over  two  GPS  stations  in  Maryland  as  it  moves  onshore.  Figure  51(a) 
shows  the  GAIM-GM  TEC  map  for  the  11-station  AFWA  grid.  Note  that  there 
are  no  other  stations  immediately  to  the  south  which  can  help  resolve  the  southern 
edge  of  the  SED.  Consequently  the  Kalman  filter  continues  to  propagate  enhanced 
densities  detected  earlier  by  the  Bermuda  GPS  station,  while  correlating  the  SED 
values  latitudinally  to  the  south.  This  results  in  the  SED  being  specified  grossly 
out  of  proportion.  However,  if  a  station  is  added  along  the  lower  East  Coast  (see 
Figure  51(b)),  GAIM-GM  detects  the  true  values  that  should  exist  over  that  area 
and  significantly  firms  up  the  solution,  reducing  the  false  TEC  values  by  over  20%. 

These  results  in  Figure  51  powerfully  illustrate  how  adding  just  one  station  to  fill 
the  gap  in  a  data-sparse  region  can  significantly  improve  GAIM-GM  TEC  results. 
They  additionally  reveal  the  disadvantage  that  this  synthetic  SED’s  trajectory  levies 
on  the  11-station  AFWA  grid  and  show  how  the  extra  station  over  Georgia  creates  a 
fairer  playing  held  for  the  numerous  cases  conducted  during  this  project. 
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Figure  51.  Differences  in  GAIM-GM  TEC  specifications  due  to  the  addition  of  a 
12th  GPS  ground  station  over  Georgia,  (a)  illustrates  how  station  sparsity  can  cause 
GAIM-GM  to  spread  perturbations  into  regions  well  beyond  the  SED.  (b)  shows  the 
correction  provided  by  the  extra  ground  station.  The  difference  between  the  two  TEC 
specifications,  shown  in  (c),  is  about  20%  over  the  new  station 
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Appendix  C.  Slant  TEC  Assimilation  Per  Time  Step  in 

GAIM-GM 


A  GPS  slant  TEC  measurement  consist  of  the  integrated  electron  densities  along 
the  path  between  a  GPS  satellite  and  a  ground  receiver.  Because  each  GPS  satellite 
is  continuously  emitting  its  electromagnetic  signal,  it  can  reach  hundreds  of  GPS 
ground  receivers  at  a  time  as  it  orbits  the  Earth.  Figure  52  illustrates  the  signal  locks 
between  GPS  satellites  and  the  AFWA  grid  at  1945  UT,  where  each  line  connecting  a 
satellite  (triangle)  to  a  ground  stations  (circle)  represents  the  signal  path  along  which 
a  slant  TEC  measurement  is  taken  and  later  ingested  by  GAIM-GM. 
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Figure  52.  Slant  TEC  paths  through  IFM  background  to  AFWA  ground  stations  at 
1945  UT.  Triangles  represent  GPS  satellites  and  circles  are  the  ground  stations.  The 
exact  number  of  slant  TEC  measurements  for  this  time  step  was  300  (not  all  slant  paths 
are  shown).  The  SED  can  be  seen  over  the  U.  S.  East  Coast 
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Table  20  lists  the  number  of  data  ingested  by  GAIM-GM  at  each  15-minute  time 
step  during  the  storm  period,  according  to  the  grid  type  and  number  of  stations  used. 
This  information  is  available  in  all  GAIM-GM  NetCDF  output  hies.  The  number 
of  slant  TEC  observations  increased  linearly  as  station  number  increased,  with  two 
exceptions.  First,  the  amount  of  data  assimilated  for  the  AFWA+400  station  grid 
were  omitted  because  the  actual  amounts  listed  in  the  NetCDF  hies  were  incorrectly 
listed  in  the  60-70  range.  The  results  in  Chapter  V  show  there  were  clearly  many 
more  data  assimilated  per  time  step.  The  second  exception  is  that  fewer  data  were 
assimilated  from  the  Ideal  720  grid  than  from  the  Ideal  360  grid.  This  suggests  that 
GAIM-GM  might  reach  data  saturation  even  before  observations  are  available  for 
every  model  grid  square,  but  future  testing  is  needed  to  verify  this  assumption. 


Table  20.  Slant  TEC  measurements  ingested  by  GAIM-GM  per  GPS  ground  station 
distribution.  Note  that  slant  TEC  counts  for  the  AFWA+400  distribution  were  not 
available  in  the  GAIM-GM  output.  This  table  pertains  to  the  modified  GAIM-GM 
version  that  trusts  data  more 


CORS  grids 

Ideal  grids 

Time  (UT) 

AFWA 

+30 

+100 

+200 

+400 

21 

90 

360 

720 

1800 

73 

284 

766 

1458 

- 

130 

549 

2201 

2031 

1815 

71 

272 

743 

1417 

- 

125 

524 

2122 

1984 

1830 

66 

254 

710 

1370 

- 

119 

499 

2018 

1798 

1845 

71 

265 

723 

1397 

- 

119 

510 

2049 

1860 

1900 

74 

279 

761 

1469 

- 

128 

538 

2140 

2042 

1915 

72 

266 

731 

1411 

- 

127 

525 

2117 

1972 

1930 

70 

252 

689 

1318 

- 

115 

497 

1993 

1644 

1945 

72 

261 

705 

1342 

- 

123 

517 

2054 

1733 

2000 

71 

273 

747 

1422 

- 

130 

535 

2125 

1865 

2015 

76 

282 

761 

1446 

- 

129 

544 

2173 

1905 

2030 

77 

273 

737 

1398 

- 

129 

533 

2140 

1908 

2045 

72 

260 

699 

1329 

- 

123 

512 

2071 

1869 
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